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ABSTRACT 


This  document  describes  a  quantitative  flow  visualization  technique  for 
transient  fluid  flow.  The  flow  visualization  technique  uses  particle  image 
velocimetry  to  obtain  instantaneous  global  measurements  of  the  fluid  flow  field. 
The  particles  in  the  flow  field  are  illuminated  with  a  multiple  pulsed  laser  sheet, 
and  the  image  of  the  illuminated  particles  is  then  captured  on  film.  The  film  is 
digitized  on  a  high  resolution  scanner,  and  the  digital  data  are  sent  to  a  computer 
system  where  either  an  autocorrelation  code  or  a  two-dimensional  fast  Fourier 
transform  code  is  sequentially  applied  to  small  regions  of  the  digitized  picture  to 
compute  the  average  displacement  of  the  particles  in  each  region.  Knowing  the 
time  between  laser  pulses  allows  the  velocity  vector  for  each  region  to  be  calculated 
from  the  displacement  data. 
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INTRODUCTION 


Past  experimental  methods  of  examining  fluid  flow  were  limited  to  point  mea¬ 
surement  techniques  such  as  laser  Doppler  velocimetry  (LDV),  hot  wire  anemom- 
etry,  and  pitot  probes.  These  techniques  were  accurate  measurement  tools  but  were 
not  useful  in  measuring  global  flow  fields.  Other  methods,  such  as  dye  injection  and 
yam  tufts,  helped  visualize  the  global  flow  field  but  did  not  provide  quantitative 
data  on  the  flow  itself.  A  tool  for  the  quantitative  measurement  of  the  global  flow 
field  needed  to  be  developed  for  hydrodynamic  investigators.  The  tool  would  have  to 
be  useful  for  steady  state  and  transient  flows. 

This  document  describes  a  flow  visualization  technique  for  transient  fluid  flow. 
Particle  image  velocimetry  (PIV)  is  used  to  obtain  instantaneous  global  measure¬ 
ments  of  the  fluid  flow  field.  The  particles  in  the  flow  field  are  illuminated  with  a 
multiple  pulsed  laser  sheet,  and  the  image  of  the  illuminated  particles  is  then  cap¬ 
tured  on  film.  The  film  is  digitized  on  a  high  resolution  scanner,  and  the  digital 
data  are  sent  to  a  computer  system  where  either  an  autocorrelation  code  or  a  two- 
dimensional  (2D)  fast  Fourier  transform  (FFT)  code  is  sequentially  applied  to  small 
regions  of  the  digitized  picture  to  compute  the  average  displacement  of  the  particles 
in  each  region.  Knowing  the  time  between  laser  pulses  allows  the  velocity  vector  for 
each  region  to  be  calculated  from  the  displacement  data. 

This  document  details  each  phase  of  the  PIV  process,  from  the  selection  of 
particles  to  the  display  of  the  computed  vectors. 


PARTICLE  SELECTION 

The  particles  in  the  flow  field  are  critical  to  the  entire  process.  The  particles 
should  be  approximately  neutrally  buoyant,  from  0.95  to  1.10  grams  per  cubic  cm, 
and  between  10  and  40  (i  in  diameter  (Katz  and  Huang,  reference  1).  Large  par¬ 
ticles  will  not  properly  track  the  flow  field  and  in  transient  flow  fields,  large  par¬ 
ticles  will  lag  the  flow  changes.  If  the  particle  is  too  small  (less  than  about  10  |i), 
the  film  will  not  be  able  to  record  its  image.  Note  that  the  image  of  the  illuminated 
particle  will  be  a  bit  larger  than  the  actual  particle  diameter  due  to  blooming  from 
the  reflected  laser  light.  Fluorescent  or  reflective  particles  can  be  used  to  enhance 
the  visibility  of  the  particles.  Fluorescent  particles  can  be  designed  to  fluoresce  at  a 
different  wavelength  than  the  laser  light.  This  will  allow  use  of  a  color  filter  to 
mask  any  unwanted  reflections  and  stray  laser  light.  Particle  reflectance  can  be 
enhanced  by  a  metal  coating  or  a  high  index  of  refraction  change  between  the  par¬ 
ticle  and  the  water. 

This  project  examined  many  particles  of  various  sizes  and  shapes,  table  1.  The 
3M  large  glass  bubbles  worked  well,  but  the  small  bubbles  were  too  small  for  the  3  x 
4  inch  area  of  interest.  The  glass  bubbles  had  to  be  washed  first  to  remove  the  silica 
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coating.  The  solid  polystyrene  particles  had  poor  reflectance  and  were  hard  to  pick 
up  on  film.  The  3  p  fluorescent  particles  were  bright  but  too  small;  they  looked 
more  like  a  dye.  The  silver  coated  particles  had  the  best  image  brightness;  however, 
at  2.6  g/cc  they  would  be  poor  at  following  the  flow.  Fluorescent  particles  of  30  p 
diameter  were  ordered  but  did  not  arrive  in  time  for  evaluation.  At  this  time  the 
best  particles  to  use  are  the  35  p  glass  bubbles. 


Table  1.  Particles  Listing 


Particle  Type 

Material 

Company 

Mean 

Diameter 

(p) 

SG 

(g/cc) 

Glass  Bubbles 

Silicon  Glass 

3M 

8.00 

1.10 

Glass  Bubbles 

Silicon  Glass 

3M 

35.00 

0.60 

Polymer  Microspheres 

Polystyrene  DVB 

Duke  Scientific 

0.5  -  24.0 

1.05 

Fluorescent  Microspheres 

Dyed  Polystyrene 

Duke  Scientific 

3 

1.05 

Metallic  Coated 

Silver  Coated  Glass 

TSI 

12 

2.6 

Plastic  Particles 

Dyed  plastic 

Dr.  J.  Katz 

35 

1.0 

SEEDING  DENSITY 

The  seeding  density  required  for  a  particular  test  is  dependent  on  size  of  the 
window  used  for  each  velocity  vector.  The  window  size  is  a  function  of  the  grid 
density  of  velocity  vectors  per  picture.  Katz  and  Huang  (reference  1)  and  Adrian 
(reference  2)  have  shown  that  both  the  FFT  and  the  correlation  methods  need  from 
6  to  10  particle  pairs  per  window  to  yield  an  accurate  velocity  vector.  Given  the  grid 
size  and  the  thickness  of  the  light  sheet,  the  required  particle  density  can  be  com¬ 
puted.  It  can  be  seen  that  for  a  given  vector  density,  as  the  viewing  area  gets 
smaller,  the  particle  density  increases  dramatically.  For  small  viewing  areas  a 
large  number  of  particles  will  cloud  the  water  and  obscure  the  light  sheet.  A  com¬ 
promise  must  be  made  between  seeding  density  and  water  clarity. 

Example: 

Window  Size: 

Experimental  Field  of  View: 

Vector  Grid: 

Thickness  of  Light  Sheet: 

Window  Volume: 

Number  of  Particles: 


0.12  cm  x  0.12  cm 

10.8  cm  x  7.2  cm,  (4.25"  x  2.83") 

90x60 
0.2  cm 

0,12  cm  x  0.12  cm  x  0.2  cm  =  0.00288  ccw 
8  particles  per  window  volume. 
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Particle  distribution  =  particles  per  window  /  water  volume  of  window 
Particle  distribution:  2777  particles/ccw. 


Particle  volume,  35  |i: 


3. 14  x  Dd  3. 14  x  0. 0035 


,-8 


2.243x  10  ccp 


Volume  of  particles  =  Particle  distribution  in  water  x  Particle  volume 
Volume  of  particles:  6.228  x  10'6  ccp/ccw. 

Density  of  particles:  1.05  g/ccp 

Particle-Water  Weight  Ratio  =  Volume  of  particles  in  water  x  Density  of  particles 
Particle-Water  Weight  Ratio:  6.540  x  10'6  g/ccw. 

Amount  of  water  in  experiment:  56.775  m3  (15,000  gal). 

Estimated  amount  of  particles  required  for  experiment : 

6.540  x  10‘5  g/ccw  x  56.775  m3  =  3.7  kg  of  particles  (8.2  lb). 

This  example  gives  a  good  estimate  for  the  amount  of  particles  required.  The 
actual  amount  will  vary  depending  on  settling,  water  clarity  and  light  sheet  design 
of  the  experiment. 


CAMERA  REQUIREMENTS 

To  get  a  respectable  viewing  area  at  high  resolution,  a  35  mm  or  larger  film 
format  must  be  used.  The  35  mm  format  has  the  advantage  over  other  formats  in 
cost  and  availability.  A  standard  35  mm  camera  can  be  used  for  single  frame  test 
shots.  A  good  quality  35  mm  camera  with  a  motorized  film  drive  will  give  from  3  to 
6  frames  per  second  (f/sec).  A  35  mm  Hultcher  camera  will  give  a  25  or  40  f/sec 
speed.  A  16  mm  LOCAM  high  speed  camera  will  give  high  frame  rates,  but  its 
viewing  area  will  be  extremely  small.  Also,  availability  of  film  in  the  16  mm  format 
is  somewhat  limited. 

Camera  shutter  speed  settings  are  a  bit  confusing.  The  setting  is  not  the  time 
required  for  the  shutter  to  open  or  close,  nor  is  it  the  time  it  takes  the  shutter  to 
traverse  the  film.  It  is  the  total  amount  of  time  for  which  any  one  section  of  the  film 
is  exposed  to  light.  For  all  standard  35  mm  cameras,  exposure  time  is  accomplished 
by  a  traveling  shutter  opening  that  passes  a  slit  of  light  across  the  film.  The  slit 
may  travel  in  the  horizontal  or  vertical  direction.  The  width  and  speed  of  the  trav¬ 
eling  slit  determine  the  exposure  time.  The  width  of  the  slit  and  speed  of  the  slit 
vary  from  camera  to  camera.  For  the  purpose  of  the  flow  visualization  project,  the 
problem  is  that  the  entire  piece  of  film  is  not  exposed  at  once  when  using  the  cam¬ 
era  speed  settings.  One  side  of  the  film  is  exposed  before  the  other. 
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For  the  flow  visualization  work,  we  want  the  picture  of  the  entire  flow  field  at 
the  same  instant  in  time.  That  means  that  the  picture  needs  to  be  taken  with  the 
shutter  in  the  fully  opened  position.  With  the  shutter  fully  opened  the  laser  system 
can  be  pulsed  to  expose  the  entire  frame  of  the  film.  To  do  this  a  strobe  synchroni¬ 
zation  signal  must  be  sent  to  the  laser  when  the  camera  shutter  is  fully  open.  On 
most  cameras  this  is  the  "X"  or  flash  setting.  A  flash  synchronization  setting  on  a 
35  mm  camera  provides  a  closed  circuit  for  the  flash  when  the  shutter  is  fully  open. 
By  connecting  a  battery  to  this  circuit  a  high  level  pulse  can  be  provided  for  other 
electronic  setups. 

The  Hultcher  camera  used  in  the  experiment  had  no  shutter  synchronization 
and  had  to  be  modified  to  give  a  pulse  before  the  shutter  opens  (figure  1).  A  slotted 
disk  was  mounted  on  the  shutter  shaft.  The  disk  has  a  slit  that  is  detected  by  an 
infrared  (IR)  emitter  detector  set.  When  the  slit  passes  the  detector  it  will  put  out  a 
5  V  transmitter  transmitter  logic  (TTL)  signal.  Since  the  slit  is  very  small  it  will  be 
a  5  V  TTL  pulse.  An  electronics  diagram  for  the  emitter  and  detector  pair  is  pre¬ 
sented  in  figure  2.  The  time  between  the  pulse  disk  output  and  the  camera  shutter 
at  the  fully  open  position  must  be  set  to  allow  ample  time  for  the  laser  system  to  get 
ready  before  the  camera  shutter  is  at  the  fully  open  position,  roughly  10  msec.  The 
pulse  delay  time  can  be  checked  by  sending  the  laser  beam,  at  low  power,  through 
the  camera  shutter  while  the  camera  is  operating.  The  laser  pulse  was  adjusted  to 
coincide  with  the  shutter  at  the  fully  open  position.  The  pulse  delay  time  was  ad¬ 
justed  as  needed  by  turning  the  pulse  disk  a  few  degrees  at  a  time. 


IR  EMITTER  IR  DETECTOR 


Figure  1.  Camera  Synchronization  Controller 
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CAMERA  DISK 


Figure  2.  Camera  Synchronization  Electronics 

For  most  tests,  the  camera  has  to  be  kept  out  of  the  flow  field.  To  get  the  re¬ 
quired  image  magnification  a  zoom  lens  was  used.  The  camera  was  kept  as  close  as 
possible  to  the  image  area.  If  the  camera  is  too  close  for  the  zoom  lens  alone,  a  lens 
and  bellows  setup  can  be  used.  A  bellows  unit  moves  the  image  plane  away  from 
the  film  plane  and  allows  a  normal  zoom  lens  to  be  used  as  a  macro  zoom  lens.  One 
must  be  careful  of  the  extremely  small  depth  of  field  when  using  the  zoom  lens  and 
bellows  arrangement. 


FILM  REQUIREMENTS 

The  viewing  area  for  the  flow  visualization  also  depends  upon  the  film  resolu¬ 
tion.  Since  color  films  typically  have  a  low  ASA  rating  and  have  resolution  lower 
than  80  lines  per  mm  (1/mm)  black  and  white  films  are  the  best  choice.  Table  2  is  a 
partial  listing  of  ASA  rating  and  resolving  powers  for  various  films  (references  3 
and  4).  While  a  25  W  laser  may  seem  powerful,  the  amount  of  light  going  to  the 
film  is  quite  small.  The  laser  pulse  is  about  3  mJ  at  8  kHz,  but  each  pulse  is  only  20 
nsec  in  duration.  So  the  light  exposure  to  the  film  is  really  quite  small.  Hence,  an 
extremely  fast  film  is  required.  Since  we  want  very  small  particles  to  properly  track 
the  flow  field,  we  need  high  resolution  film  to  capture  the  particles  accurately.  Tri- 
X  Pan  with  an  ASA  rating  of  400  was  tested  and  was  found  to  be  too  slow.  Based  on 
this  test  and  the  capabilities  of  the  film,  KODAK  T-MAX  P3200  black  and  white 
film  was  tested.  This  film  has  an  ASA  rating  of  3200  and  a  resolution  of  125  1/mm. 
The  film  worked  well  with  the  particles  tested. 


Table  2.  Film  Listing 


Film  Type 

ASA  Rating 

Resolving  Power  (1/mm) 

BLACK  AND  WHITE 

Technical  Pan  2415 

25 

320 

Plus-X  Pan 

125 

125 

Tri_X  Pan  (TX) 

400 

100 

T-MAX  P3200 

800  -  3200 

125 

COLOR  FILMS 

Ectar  1000 

1000 

80 

Gold  Extapress 

1600 

40-80 

Due  to  the  resolution  of  the  film,  the  size  of  the  particle  must  be  taken  into 
account  when  setting  up  the  camera  field  of  view.  The  particle  should  be  a  mini¬ 
mum  of  3  to  4  lines  wide  on  the  negative.  It  is  not  easy  to  calculate  the  maximum 
field  of  view  for  a  given  particle.  While  the  particle  has  a  given  size,  its  image  on 
film  can  be  larger  due  to  the  reflective  properties  of  the  particle,  clarity  of  the  water, 
and  laser  power.  Testing  for  each  particle  would  have  to  be  conducted  to  get  its 
image  size  on  film.  In  the  experimental  setup  used  for  this  project  (2  mm  light 
sheet  with  the  laser  at  25  W  and  camera  lens  setting  of  f/5.6),  the  12  p  silver  coated 
particle  produced  a  3  line  picture  with  a  field  of  view  of  10.8  cm  x  72  cm  (4.25  inch  x 
2.83  inch).  For  the  same  setup  the  35  p  glass  bubbles  produced  a  4  line  picture  of 
the  particle. 

One  other  particle  parameter  is  the  brightness  of  the  image.  The  particle 
needs  to  produce  a  high  contrast  image  in  order  for  the  computer  to  "see"  the  par¬ 
ticle. 


LASER  SELECTION 

The  PIV  technique  requires  a  high  power,  high  speed,  pulsed  laser  source. 
Laser  options  are  limited  to  a  continuous  wave  (CW)  laser  with  a  chopping  wheel, 
(Nd:Yag)  or  a  metal  vapor  laser.  The  selected  laser  must  be  able  to  output  short 
pulses,  10  nsec,  at  high  repetition  rates,  2  to  10  kHz.  The  CW  laser  with  a  chopper 
wheel  is  limited  to  4  kHz  with  a  relatively  large  pulse  width  of  125  microseconds. 
The  Nd: Yag  laser  runs  around  20  to  200  pps  with  a  pulse  width  of  7  nsec  with  per 
pulse  power  of  around  120  mJ.  While  the  Nd:Yag  provides  a  high  power  pulse,  its 
use  requires  two  lasers  to  give  double  pulsed  capability.  For  this  project  the  best 
compromise  between  high  repetition  rate  and  power  per  pulse  was  a  25  W  copper 
vapor  laser.  This  laser  uses  copper  metal  vaporized  at  extremely  high  temperatures 
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and  voltages.  The  copper  vapor  lases  and  emits  a  green  yellow  beam.  The  25  W  CJ 
laser  used  in  this  testing  puts  out  25  W  average  power  at  a  nominal  8  kHz  with  a 
peak  power  of  3.5  mJ  at  8  kHz  and  2  mJ  at  2  kHz.  The  laser  pulse  rate  range  is 
from  2  to  10  kHz. 

To  keep  the  copper  metal  vaporized  and  lasing,  the  laser  tube  must  be  kept 
hot.  To  keep  the  laser  hot,  it  must  be  running  near  its  design  frequency,  in  this  case 
8  kHz.  For  the  PIV  work  laser  pulse  frequency  may  be  as  low  as  2  kHz.  However, 
due  to  the  starting  and  stopping  required  by  the  camera,  20  or  40  Hz,  the  laser  off 
time  would  allow  the  system  to  cool  down.  A  laser  interface  box,  ESC-2000,  which 
will  run  the  laser  at  8  kHz  normally  and  will  keep  the  laser  running  at  the  required 
8  kHz,  stop  it,  run  it  at  the  requested  frequency  and  then  start  it  again,  is  connected 
to  the  laser.  A  laser  start  pulse  will  tell  the  laser  to  stop  running  at  8  kHz,  open  the 
laser  shutter,  EMS-2000,  and  run  at  the  external  trigger  frequency.  This  takes 
about  5  msec.  A  stop  pulse  command  will  tell  the  laser  to  stop  running  at  the  exter¬ 
nal  trigger  frequency,  close  the  shutter  and  run  at  the  8  kHz  again.  This  takes 
about  5  msec. 


LIGHT  SHEET  DESIGN 

PIV  works  by  taking  a  picture  of  particles  in  the  flow  field.  The  selection  of 
which  particles  to  view  is  made  by  forming  the  laser  light  into  a  light  sheet  and 
projecting  the  light  sheet  into  the  water  (figure  3).  The  camera  is  held  perpendicu¬ 
lar  to  the  light  sheet  for  viewing  the  flow  field.  The  light  sheet  must  be  thin  in 
order  to  reduce  errors  caused  by  particles  going  at  angles  to  the  investigated  path. 

Due  to  the  large  output  beam  diameter  (20  mm)  of  the  copper  vapor  laser,  the 
laser  beam  must  be  focused  to  a  smaller  size  (about  2.5  mm).  The  small  beam  is 
then  passed  through  a  cylindrical  lens  to  form  a  light  sheet.  All  of  the  lenses  should 
be  set  to  minimize  reflections  back  into  the  unstable  resonator.  Reflections  into  the 
unstable  resonator  will  contaminate  the  quality  of  the  generated  light. 

For  this  testing  the  beam  was  focused  by  passing  it  through  a  2  inch  f/1000.0 
mm  lens.  An  f/125.0  mm  lens  was  then  placed  after  the  focusing  point  to  refocus 
and  collimate  the  smaller  beam.  The  small  beam  was  then  passed  through  the 
cylindrical  lens  to  form  the  light  sheet. 
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PLANO-CONVEX  LENS 
NEWPORT  LENS  KPX217 
F  1000.0  mm 


20  mm  BEAM 


MIRROR 


20  mm  BEAM 


LIGHT  SHEET 


PLANO-CONVEX  LENS 
NEWPORT  LENS  KPX190 
F  125.0  mm 


CYLINDRICAL 
PLANO-CONVEX  LENS 
NEWPORT  LENS  CKX025 
F  25.40  mm,  R  13.13  mm 


AA-2000  ATTENUATOR 

EMS-2000  ELECTROMAGNETIC  SHUTTER 

BD90-2000  BEAM  STEERING  OPTICS 


Figure  3.  Optical  Train 


LASER  PULSE  CONTROL 

The  laser  pulse  must  be  controlled  so  that  the  laser  will  pulse  when  the  camera 
is  at  the  fully  opened  position.  The  camera  sends  out  a  ready  pulse  for  the  com¬ 
puter  to  send  a  signal  to  pulse  the  laser.  A  Masscomp  computer  was  used  to  control 
the  pulsing  of  the  laser.  The  computer  front  panel  clocks  (figure  4)  and  the  pro¬ 
gram  laser_control.c  are  used  to  control  the  laser  start  and  stop  pulses  and  the 
external  trigger  pulses.  The  computer  program  sets  up  the  delay  time  between  the 
camera  ready  synchronization  and  the  laser  start  synchronization.  The  program 
must  be  told  the  proper  time  settings  for  the  shutter  delay  line,  PIV  frequency,  and 
number  of  pulses.  The  code  is  currently  set  to  run  the  laser  at  one  frequency.  By 
using  the  transient  subroutine  in  the  code,  the  clock  can  vary  the  laser  frequency 
during  the  test.  To  eliminate  false  triggering  of  the  laser  due  to  electronic  noise,  the 
program  sets  the  trigger  level  on  the  Masscomp  clock  to  2.5  V.  This  eliminates 
noise  problems  but  does  cause  some  problems  with  other  Masscomp  users.  Once  the 
program  is  finished,  the  clocks  are  reset  to  the  standard  1 V  setting.  The  program 
also  puts  out  a  graphic  file  that  displays  the  clock  pulse  setup  and  timing.  An  ex¬ 
ample  of  the  timing  chart  is  shown  in  figure  5. 
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LASER  START  GATE  OUTPUT  PULSE 

3£) 

CAMERA  READY  INPUT  PULSE 

C0\  O  ~^r 

LASER  STOP  GATE  OUTPUT  PULSE 

_l£>  O 

LASER  PIV  BURST  OUTPUT 

c0C 

*© 

6W)  O^r 

cQr  O-^r 

Figure  4.  Masscomp  CK10  Front  Panel 


Figure  5.  Masscomp  Clock  Timing  Chart 
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CLOCK  SYNCHRONIZATION  OPERATIONS 


The  camera  ready  pulse  will  come  from  the  camera  ready  line  and  go  into  the 
CK3,  4  and  5  C  inputs  of  the  CK10.  The  camera  ready  pulse  is  prior  to  the  camera 
shutter  fully  open  position  and  is  set  to  account  for  the  camera  disk  offset,  allow 
time  for  the  laser  to  stop  pulsing  at  8  kHz,  and  open  the  laser  shutter,  about  5  msec. 
The  CK3  output,  start  laser  pulse,  will  tell  the  laser  to  stop  lasing,  open  the  shutter, 
and  pulse  at  the  external  trigger  frequency.  The  start  laser  pulse  is  delayed  to 
allow  time  for  the  laser  shutter  to  open  and  at  the  same  time  it  must  be  set  to  be 
coincident  with  the  camera  shutter  fully  open  position.  This  is  done  by  trial  and 
error.  The  CK4  output  is  set  to  send  a  close  laser  pulse  once  the  laser  has  finished 
its  pulsing.  The  CK5  output  will  tell  the  computer  when  to  pulse  and  for  how  long; 
it  is  connected  into  the  CK6  C  input.  The  CK  6  output  is  connected  to  the  laser 
external  trigger  input.  The  external  trigger  frequency  is  determined  by  the  ffbw 
velocity  and  viewing  area  of  the  experiment. 


AUTOMATIC  SCANNING  SYSTEM 

An  automatic  film  scanning  system  was  designed  around  a  Nikon  scanner. 
The  system  is  controlled  by  software  on  the  PC,  which  controls  the  scanning,  pro¬ 
cessing,  and  film  movement  (figure  6).  The  scanning  is  performed  on  a  Nikon  Ls- 
3510  AF  high  resolution  scanner.  The  film  is  moved  by  two  Newport  Model  850 
linear  actuators  and  one  Model  495  rotary  actuator.  The  film  movement  is  con¬ 
trolled  by  the  PC  via  a  Newport  PMC-300  Power/Interface  Module. 


Figure  6.  Film  Transport  System 


Prior  to  running  the  auto  scan  software,  one  must  align  the  film  on  the  film 
holder  and  the  film  transport  pins  must  be  aligned  with  the  film  sprocket  holes. 
The  film  transport  alignment  is  accomplished  by  using  the  Window  for  Actuators 
and  Rotary  stages  Program  (WARP)  code  provided  with  the  PMC-300  control  sys¬ 
tem.  It  must  also  be  verified  that  the  data  analysis  software  is  compiled  and  run¬ 
ning  on  the  Cray  and  that  the  Cray  is  on-line. 


CONTROL  CODE 

The  runpiv.c  program  controls  the  digitizing  and  submission  process  by  call¬ 
ing  other  computer  codes.  The  runpiv.c  code  will  ask  for  a  run  id,  starting  frame, 
and  number  of  frames.  It  will  use  this  information  to  sequence  through  the  film 
and  create  appropriate  data  files.  The  files  will  use  the  format  of 
idrun_numberfrframe_number.file_type,  where  file  type  is  .job  for  the  Cray  job 
file,  .snd  for  the  Cray  transfer  shell  file,  .put  for  the  PC  file  transfer  file,  and  .dsp 
for  the  computed  displacement  data  file. 


FILM  SCANNING 

The  control  code  calls  the  scanpic.c  program  to  scan  the  image  from  the  film. 
The  image  is  scanned  on  a  NIKON  LS-3510AF  color  scanner  with  a  4500  x  3000 
pixel  resolution.  The  image  is  scanned  in  as  256  gray  shades.  Control  of  the  scan¬ 
ner  is  by  the  IEEE-488  bus  from  the  PC.  The  resulting  file  is  stored  with  a  .img 
suffix.  The  image  file  is  rather  large  (27  Mb). 


JOB  SUBMISSION 

The  control  code  next  runs  submit.bat  to  transmit  the  image  and  the  analysis 
job  to  the  Cray.  The  transmission  is  made  using  Ethernet  TCP/IP  protocol.  The 
user  must  have  a  current  account  on  the  Cray  computer  and  have  the  ftp  network 
files  set  up  for  this  account.  Transfer  of  the  image  file  to  the  Cray  can  take  as  long 
as  1  hour.  Files  are  automatically  removed  from  the  Cray  after  the  job  is  completed. 
After  data  analysis  is  completed  the  displacement  data  file  is  sent  to  the  Silicon 
Graphics  Incorporated  workstation  for  post-processing  and  display. 


FILM  ADVANCE 

The  control  code  then  calls  advance.c  to  move  the  film  transport  to  advance 
the  film  to  the  next  frame.  The  transport  comprises  two  linear  actuators  mounted 
to  engage  and  move  the  film  with  two  pins  (figure  6).  A  rotary  actuator  is  also 
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mounted  to  wind  the  film  during  the  process.  The  transport  pins  are  moved  into  the 
holes  on  the  film  and  then  the  film  is  pulled  one  frame  length.  The  pins  are  pulled 
away  from  the  film  and  the  transport  is  moved  back  to  the  starting  position. 


PARTICLE  IMAGE  VELOCIMETRY 

Particle  image  velocimetry  uses  multiple  exposures  of  moving  particles  to 
measure  particle  velocity.  By  using  a  pulsed  laser  to  illuminate  the  particle  at 
different  times  and  recording  the  particle  images,  one  can  measure  the  distance 
between  the  particles.  Knowing  the  time  between  pulses  allows  the  velocity  to  he 
computed.  To  visualize  this,  consider  a  laser  pulse  at  time  t0.  Images  of  the  par¬ 
ticles  are  recorded  on  film.  The  laser  is  pulsed  again  some  time  later,  time  tx;  the 
particles  have  moved  a  distance  d,  which  depends  on  the  particle  velocity  v  and  the 
time  between  pulses  dt.  The  image  of  the  particles  at  time  tx  will  be  shifted  from 
the  image  of  the  particles  at  time  t0.  Both  exposures  are  recorded  on  one  piece  of 
film  (figure  7). 


Figure  7.  Double  Exposed  Particle  Image 

The  particle  velocity  can  be  computed  directly  from  the  picture.  A  large  picture 
with  many  particles  is  broken  down  into  small  windows.  Each  window  must  be 
sized  so  that  all  of  the  pulsed  images  from  a  particle  fit  in  that  window.  At  large 
flow  speeds  this  can  require  large  windows.  The  laser  frequency  can  be  increased  to 
reduce  the  window  size,  but  this  does  increase  the  minimum  measurable  particle 
velocity.  The  average  velocity  of  each  window  is  computed  and  displayed  by  the 
computer. 

There  are  two  primary  methods  for  computing  the  average  velocity  of  the 
particles.  The  first  method  is  a  spatial  autocorrelation  method  that  uses  overlaying 
maps  of  the  image  and  shifts  the  maps  to  determine  the  best  correlation.  The  sec- 


ond  method  uses  two  2D  FFT’s  to  compute  the  spatial  frequency  of  the  particle 
pairs. 

The  image  is  broken  into  small  windows  to  compute  the  average  particle  veloc¬ 
ity  for  each  window.  The  size  of  the  window  must  be  larger  than  the  maximum 
expected  particle  displacement.  For  high  accuracy  and  wide  range  the  window 
should  be  as  large  as  possible.  Unfortunately,  large  windows  will  lower  the  number 
of  velocity  vectors  and  grid  resolution. 


AUTOCORRELATION  METHOD 

The  autocorrelation  method  uses  image  shifting  to  find  the  best  correlation. 
This  is  accomplished  by  taking  a  copy  of  the  double  exposed  image  and  placing  that 
copy  directly  over  the  original  image  (Katz  and  Huang,  reference  1).  The  initial 
position  of  the  copy  is  at  the  0,0  point.  Now  consider  that  the  image  is  in  digital 
format.  For  simplicity  let  the  dark  pixels  be  1  and  the  clear  points  be  0  (the  gray 
shade  of  the  digitized  image  is  used  in  the  computer  program).  Go  through  each 
pixel  of  the  original  image,  take  its  value,  and  multiply  it  by  the  value  of  the  pixel 
on  the  copy  directly  above  it.  Do  this  for  all  of  the  pixels  and  add  the  values  to¬ 
gether.  This  will  give  the  autocorrelation  value  for  the  zero  position.  This  will  be 
the  peak  autocorrelation  value,  since  it  has  the  best  correlation  (figure  8).  If  the 
copy  on  the  top  is  moved  by  a  single  pixel  in  both  the  x  and  y  directions,  and  the 
process  is  repeated  for  each  movement,  the  new  value  will  be  the  autocorrelation 
value  for  that  pixel  displacement.  The  peak  value  from  the  zero  position  is  com¬ 
pared  with  the  next  highest  peak.  The  next  highest  peak  will  be  for  the  position 
when  the  t0  particles  on  the  copy  best  overlay  the  tx  particles  on  the  original.  The 
distance  between  the  zero  position  and  the  next  highest  peak  is  the  average  dis¬ 
placement  for  those  particles.  This  is  the  basic  process  behind  autocorrelation  PIV. 
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The  autocore.c  program  was  written  to  take  a  multiple  pulsed  exposure  of  the 
particle  image  and  perform  the  autocorrelation.  The  program  reads  in  the  digitized 
picture  and  breaks  up  the  picture  into  windows.  The  average  displacement  for  each 
window  is  computed.  Each  window  should  be  sized  to  have  about  8  particle  pairs  in 
it  for  the  autocorrelation  to  be  accurate.  A  smaller  window  will  run  faster  and  give 
more  vectors,  but  the  accuracy  of  the  vector  will  be  poorer  since  the  number  of  * 
particle  pairs  will  be  lower.  A  large  window  will  be  more  accurate  but  will  take  very 
long  to  compute  and  will  yield  few  velocity  vectors.  The  size  of  the  window  will  have 
to  be  determined  by  looking  at  the  negative.  The  program  also  has  to  know  how  big, 
in  pixels,  the  particle  image  appears.  This  forces  the  code  to  go  past  the  zero  over¬ 
lap  point,  which  is  the  diameter  of  the  particle  in  pixels. 

The  particle  displacement  vector  is  computed  by  finding  the  distance  from  the 
zero  point  to  the  location  of  the  next  highest  correlation  peak.  The  location  of  that 
peak  may  not  provide  an  accurate  measurement  of  the  average  displacement  of  the 
particles  in  the  window.  To  get  a  better  average  displacement  value,  the  average 
location  value  of  the  peak  and  its  surrounding  area  is  computed.  This  is  performed 
in  the  code  by  averaging  the  correlation  weighted  locations  to  compute  the  best 
average  peak  location. 

The  code  will  go  through  each  window  and  compute  the  autocorrelation  for  the 
positive  x  and  y,  and  for  the  positive  x  and  negative  y,  and  output  the  vectors  into  a 
data  file. 

This  code  takes  a  considerable  amount  of  time  (1  to  4  hours)  on  the  Cray  to 
compute  all  of  the  vectors.  If  one  had  sufficient  knowledge  of  the  flow  field  one 
could  reduce  the  time  of  computation  by  giving  the  code  good  initial  starting  points 
and  limited  velocity  bands. 


TWO-DIMENSIONAL  (2D)  FAST  FOURIER  TRANSFORM  (FFT)  METHOD 

The  other  method  to  compute  particle  displacements  uses  2D  FFT  (Adrian, 
reference  2).  Given  a  window  of  multiple  particle  exposures  (figure  9),  a  2D  FFT 
can  be  used  to  get  the  average  spatial  displacement  between  particles.  A  2D  FFT  is 
performed  on  a  copy  of  the  original  image;  a  display  of  the  magnitude  output  will 
show  the  interference  fringes  (figure  10),  identical  to  the  optically  produced  Young’s 
fringes  .  A  2D  FFT  on  the  interference  fringes  will  produce  a  peak  at  the  average 
displacement  point  (figures  11  and  12).  Both  the  highest  peak  and  the  second  high¬ 
est  peak  are  found.  The  ratio  of  the  highest  to  the  second  highest  peak  can  be  used 
as  an  indicator  of  the  quality  of  the  data.  A  ratio  of  1.25  or  higher  has  been  shown 
to  reflect  good  data  (Keane  and  Adrian,  reference  5). 
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Figure  9.  Simulated  Multiple  Laser  Pulsed  Particle  Field 


Figure  10.  Computed  Young's  Fringes 
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Figure  11.  Two-Dimensioned  Display  of  Peak  Correlation 


Figure  12.  Three-Dimensional  Display  of  Peak  Correlation 
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Limitations  do  exist  for  the  2D  FFT  method.  The  image  window  must  encom¬ 
pass  about  8  particle  pairs  to  be  accurate  and  the  order  of  the  FFT  must  be  greater 
than  half  the  greatest  particle  displacement,  also  known  as  the  Nyquist  criterion. 
High  order  FFTs  can  be  implemented  by  zero  padding  of  the  image  window.  There 
are  two  versions  of  the  code  using  the  2D  FFT  method.  The  first  uses  a  C  version  of 
the  FORTRAN  2D  FFT  code  in  the  Numerical  Recipes  book.  This  code  can  be  rfin 
on  the  SGI,  taking  about  2  hours,  or  on  the  Cray,  taking  about  1  hour.  Using  the 
Cray  intrinsic  2D  complex  FFT  subroutine  reduces  the  computation  time  to  around 
5  minutes. 


DIRECTIONAL  AMBIGUITY 

Since  PIV  is  represented  as  a  series  of  dots  on  the  film  there  is  no  information 
as  to  actual  particle  direction.  If  the  direction  of  the  flow  is  known  in  one  axis,  the 
direction  can  be  determined  for  the  other  axis.  If  no  direction  is  known,  as  is  the 
case  in  some  flows,  the  computer  will  be  able  to  determine  angle  but  not  the  direc¬ 
tion  of  the  flow.  Adrian  (reference  6)  has  shown  that  this  ambiguity  can  be  resolved 
by  a  technique  called  image  shifting.  An  additional  displacement,  larger  than  the 
most  negative  particle  displacement,  is  added  to  the  particles  by  moving  the  camera 
or  the  film  during  the  filming  process.  Landreth  and  Adrian  (reference  7)  developed 
a  technique  that  uses  Pokel  cells  for  electro-optical  movement  of  the  image.  The 
artificial  displacement  is  removed  during  processing  and  the  direction  can  be  re¬ 
solved.  While  this  method  provides  an  increased  vector  window  it  yields  a  smaller 
number  of  vectors  per  frame  and  a  smaller  field  of  view.  For  most  flows  a  high 
vector  density  is  required  and  the  direction  can  be  determined  from  other  sources. 


DATA  PRESENTATION 

After  the  displacement  vectors  are  computed  they  can  be  displayed  on  the  SGI 
IRIS.  The  vector.c  code  will  display  the  vectors  on  the  screen  and  also  dump  an 
HPGL  file  for  plotting  the  vectors.  The  vector  code  can  dump  all  of  the  vectors  or 
try  to  eliminate  the  extraneous  vectors.  The  code  uses  several  methods  to  eliminate 
the  bad  vectors.  The  first  method  is  to  use  the  ratio  of  the  highest  to  second  highest 
peak.  The  ratio  can  be  adjusted  to  balance  the  quality  of  the  data  with  the  amount 
of  the  data.  The  code  also  loops  through  each  computed  vector  and  compares  the 
magnitude  and  angles  of  the  selected  vector  to  those  at  the  adjacent  vectors.  If 
there  are  not  at  least  two  adjacent  vectors  within  a  given  percentage  of  the  selected 
vector,  the  selected  vector  is  not  displayed. 
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EXPERIMENTAL  SETUP 


A  small  test  chamber  was  set  up  to  test  the  PIV  system  (figure  13).  The  cham¬ 
ber  was  kept  small  to  reduce  the  amount  of  particles  used  per  test.  A  10-gallon 
glass  aquarium  tank  was  used  to  hold  water.  A  small  foot-operated  bellows  pufnp 
pushed  water  out  of  a  5/8  inch  copper  tube.  The  light  sheet  was  approximately  0.2 
cm  wide  and  was  aligned  with  the  tube  exit.  A  35  mm  Nikon  F2  camera  was  used 
with  a  200  mm  zoom  lens  on  a  bellows  extension.  The  camera  was  set  approxi¬ 
mately  4  feet  from  the  copper  tube.  The  field  of  view  was  108  cm  x  72  cm. 


WATER  TANK 


Figure  13.  Experimental  Setup 


18  . 


The  laser  was  always  run  at  maximum  power  (25  W).  The  typical  laser  fre¬ 
quency  was  2  kHz.  Various  f  stop  settings  were  tried.  Low  settings  produced  a  good 
image  intensity  but  a  poor  depth  of  field;  high  settings  limited  the  amount  of  light  to 
the  film.  The  best  setting  for  this  setup  was  f/5.6. 

A  number  of  experiments  were  tried  to  find  the  best  number  of  pulses  per 
frame.  The  minimum  required  for  PIV  is  2.  Three  pulses  typically  gave  better' 
results.  When  using  more  than  3  pulses  per  frame  the  image  tended  to  loose  clarity 
due  to  the  increased  background  image  noise.  A  setting  of  3  pulses  per  frame  is 
recommended  for  PIV. 

The  Hultcher  camera  was  tested  at  25  f/sec.  The  film  was  not  put  through  the 
beater  spool  so  that  the  film  would  be  constantly  moving.  This  film  movement 
provided  image  shifting  for  resolving  the  directionality  of  the  flow.  The  window  size 
must  be  enlarged  to  fit  the  increased  displacement  of  the  particles.  Unfortunately, 
the  larger  window  size  yields  fewer  vectors.  The  camera  was  run  in  its  normal 
mode  to  record  the  jet  flow  at  25  f/sec.  Unfortunately,  the  film  was  not  developed 
properly  and  there  was  no  time  left  in  the  program  to  repeat  the  test. 


EXPERIMENTAL  RESULTS 

The  PIV  system  provided  global  velocity  information  on  the  flow  field.  The 
first  test  (figure  14)  shows  a  jet  of  water  starting  up  in  a  quiescent  tank;  flow  direc¬ 
tion  is  in  the  positive  y  direction.  All  data  are  displayed  with  no  effort  taken  to 
remove  the  extraneous  (bad)  vectors.  These  data  were  collected  using  a  low  seeding 
density  of  the  12  p.  silver  coated  particles.  The  laser  was  set  at  10  kHz  producing  25 
W.  Ten  pulses  were  used  with  a  2  mm  thick  light  sheet.  Figure  15  shows  the  data 
after  the  image  has  been  cleaned  up.  The  water  plume  and  vortices  generated  at 
startup  are  evident  in  both  graphs. 

The  second  test  shows  a  starting  flow  with  3  pulses  at  2  kHz  (figures  16  and 
17).  The  startup  vortices  are  evident  in  figure  18.  In  figure  16  the  indication  of 
vortices  is  lost  due  to  the  removal  of  the  particle  diameter  limited  displacement 
vectors.  The  magnification  of  the  image  onto  the  film  must  be  sized  to  make  sure 
that  the  small  velocities  are  not  lost  in  the  image  noise.  The  third  test  shows  a 
steady  state  jet  flow,  with  raw  data  in  figure  18  and  reduced  data  in  figure  19. 


FILE  NAME  :  frfflAO.dsp 
RUN  ID  NUMBER  :  0 
Frese  Nuder  :  0 

Picture  size  :  4500  X  3000  pixels 
Window  size  :  40  X  40  pixels 
Vector  Count  :  112  X  75  :  total  B400 
Conversion  :  1  pixel  ■  0.052000  ft/sec 


Figure  14.  Startup  Jet  Flow,  10  kHz  Laser  Pulse  (Raw  Data) 
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Figure  15.  Startup  Jet  Flow,  10  kHz  Laser  Pulse  (Reduced  Data) 
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Figure  16.  Startup  Jet  Flow ,  2  kHz  Laser  Pulse  (Raw  Data) 
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Figure  17.  Startup  Jet  Flow,  2  kHz  Laser  Pulse  (Reduced  Data) 
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Figure  18.  Steady  State  Jet  Flow,  2  kHz  Laser  Pulse  (Raw  Data) 
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Figure  19.  Steady  State  Jet  Flow ,  2  kHz  Laser  Pulse  (Reduced  Data) 
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All  of  the  data  show  good  performance  at  the  high  velocity  regions  of  the  flow. 
Since  the  testing  was  with  low  speed  flow,  some  of  the  lower  magnitude  vectors  fall 
below  the  minimum  measurable  velocity  vector  limit  and  are  not  displayed.  The 
velocity  limit  is  an  artifact  of  the  size  of  the  particle  image  on  the  film  and  is  dis¬ 
cussed  in  the  next  section. 


ACCURACY 

The  accuracy  of  the  PIV  system  is  dependent  on  a  number  of  variables.  The 
optical  magnification  and  distortion  can  be  measured  by  photographing  a  grid  pat¬ 
tern  in  the  test  chamber.  A  good  measurement  of  magnification  is  required  to  get 
an  accurate  vector  displacement  value.  Accurate  tracking  of  the  flow  can  be  best 
achieved  by  a  neutrally  buoyant  30  to  40  p  diameter  particle.  The  low  limit  of  mea¬ 
sured  velocity  is  the  minimum  possible  measured  displacement.  Since  the  particles 
must  be  separated  by  at  least  one  pixel,  the  minimum  possible  measured  displace¬ 
ment  is  equal  to  the  particle  image  diameter  plus  one.  The  vector  window  is  sized 
by  the  maximum  velocity,  laser  pulse  frequency,  and  required  accuracy. 

Example: 

Real  Image  Area:  10.8  cm  x  7.2  cm  = 

Window  Size:  0.12  cm  x  0.12  cm 

Particle  diameter  =  3  pixels 
Pixel  ratio  =  0.0024  cm/pixel 
Number  of  laser  pulses  per  frame  =  2. 

Maximum  Flow  Velocity:  10  m/sec  =  1000  cm/sec 

Set  particle  displacement  for  the  maximum  flow  velocity  to  40  pixels  to  fit 
within  the  50  pixel  wide  window. 

Measurement  accuracy  =  1  /  particle  displacement  =  1  /  40  =  2.5%  of  span  accuracy 
For  this  case  the  measurement  accuracy  would  be  ±  12.5  cm/sec. 

Laser  Pulse  Rate  =  flow  rate  /  (pixel  ratio  x  particle  displacement) 

=  (1000  cm/sec  /  (0.0024  cm/pixel  x  40  pixels)  =  10,416  Hz. 

Increasing  the  number  of  laser  pulses  per  frame  would  require  larger  windows  or 
higher  laser  pulse  rates. 

Minimum  measurable  velocity  =  (particle  diameter  +  1)  x  maximum  flow  rate  / 
particle  displacement  =  (3  +1)  x  1000  cm/sec  /  40  =  100  cm  /sec. 

This  example  show  values  for  high  vector  density  with  reasonable  accuracy. 
Larger  windows  may  be  required  to  get  higher  accuracy.  Depending  on  the  particle 
density  in  the  water,  larger  windows  may  also  be  required  to  get  8  particle  pair 


=  4500  x  3000  pixels 
=  50  x  50  pixels 


images  in  each  window.  The  scanning  of  the  film  is  considered  to  be  highly  accu¬ 
rate.  The  computation  of  the  vector  by  the  FFT  method  can  vary  in  accuracy  de¬ 
pending  on  the  computer  system  used. 


CONCLUSIONS 

Particle  image  velocimetry  can  provide  the  required  global  velocity  field  flow 
measurement  tool  required.  Particle  seeding  requirements  must  be  adjusted  to 
each  test.  Laser  pulse  rates  depend  on  the  highest  flow  speed  and  required  vector 
density.  The  lowest  measurable  displacement  is  also  dependent  on  the  laser  pulse 
rate.  Directional  ambiguity  of  the  vector  is  a  trade-off  with  the  vector  density. 
Most  of  the  foreseen  work  requires  a  high  vector  density. 


RECOMMENDATIONS  AND  FUTURE  WORK 

While  the  current  PIV  system  works,  improvements  can  be  made.  The  auto¬ 
loading  system  could  be  improved  by  cutting  the  film  into  slides  and  using  an  auto¬ 
matic  slide  loader  on  the  scanner.  The  overall  speed  of  the  system  could  also  be 
increased  by  reducing  the  transfer  time  to  the  Cray  or  doing  the  computations 
locally.  A  dedicated  array  processor  on  the  PC  would  increased  the  efficiency  of  the 
system.  Humphreys  (reference  8)  has  shown  that  higher  particle  density  could  be 
achieved  by  using  the  computed  average  displacement  vector  for  each  window  to 
find  every  particle  pair  displacement  vector  in  the  window.  This  technique  could  be 
added  to  the  PIV  software. 


REFERENCES 


1.  J.  S.  Katz  and  T.  T.  Huang,  "Quantitative  Visualization  of  External  and  Inter¬ 
nal  Flows  by  Implementing  the  Particle  Displacement  Velocity  Method,"  Ex¬ 
perimental  and  Numerical  Flow  Visualization.  FED-vol  128,  ASME  1991. 

2.  R.  J.  Adrian,  "Multi-Point  Optical  Measurements  of  Simultaneous  Vectors  in 
Unsteady  Flow:  A  Review,"  Int.  J.  Heat  Fluid  Flow.  7,127  (1986). 

3.  "Kodak  Black  and  White  Films  Data,  Darkroom  Handbook,"  Eastman  Kodak 
Company,  Rochester,  NY,  14650. 

4.  "Kodak  T-Max  Professional  Films,"  Kodak  Publication  No.  F-32,  Eastman 
Kodak  Company,  Rochester,  NY,  14650. 

5.  R.  D.  Keane  and  R.  J.  Adrian,  "Optimization  of  Particle  Image  Velocimeters 
with  Multiple  Pulse  Systems,"  Proceedings  of  the  5th  International  Symposium 
on  Applications  of  Laser  Techniques  to  Fluid  Mechanics.  Lisbon,  Portugal, 
paper  12.4,  1990. 

6.  R.  J.  Adrian,  "Image  Shifting  Technique  to  Resolve  Directional  Ambiguity  in 
Double-Pulsed  Velocimetry,"  Applied  Optics,  vol.  25,  no.  21, 1  November  1986. 

7.  C.  C.  Landreth  and  R.  J.  Adrian,  "Electrooptical  Image  Shifting  for  Particle 
Image  Velocimetry,"  Applied  Optics,  vol.  27,  no.  20, 15  October  1988. 

8.  W.  M.  Humphreys,  "A  Histogram-Based  Technique  for  Rapid  Vector  Extraction 
from  PIV  Photographs,"  Fourth  International  Conference  on  Laser  Anemom- 
etrv  -  Advances  and  Applications.  Cleveland,  OH,  1991. 


28  . 


Appendix  A 


LASER  PULSE  CONTROL  CODE:  LASER.CON.C 


A-l/A-2 
Reverse  Blank 


/*  Laser  Control  Program  for  ESC-2000  on  the  */ 

/*  Copper  Vapor  Laser  */ 

/*  Kenneth  LaPointe  ,  NUWC,  Code  8322  */ 

/*  last  Rev  Date :  4/23/92  */ 

/*  Compile  with  cc  laser_con.c  -lmr  -lgp  -lm  -o  laser_con*/ 

/*  Camera  sync  pulse  comes  is  TTL  pluse  and  comes  from 
the  35  mm  camera.  The  sync  pulse  is  connected  to  the  “C” 
inputs  of  clocks  3,4  and  5.  The  output  of  5  is  connected  into 
the  “C”  input  of  clock  6. 

Clock  3  outputs  the  laser  start  gate. 

Clock  4  outputs  the  laser  stop  gate. 

Clock  5  outputs  a  high  level  trigger  to  clock  6. 

Clock  6  outputs  the  laser  pulse  signals  to  the  Laser 


#include  </usr/include/mr.h> 

#include  <stdio.h> 

#include  <math.h> 

#include  <libmp.h> 

long  gls[SIZEOFGCA]; 

int  ij,k,l; 

int  clock_path[10]; 
intclk_da[10]; 
int  nclocks; 

int  sixmhz  =  11;  /*  six  Mhz  rising  edge 

int  highlevelgate  =  4; 
int  lowlevelgate  =  5; 

int  risingedgegate  =  6  ;/*  Rising  Edge  Gate  ,  C  input  of  this  clock 

int  fallingedgegate  =  7; 

int  startlow  =  0; 

int  starthigh  =  1; 

int  clockmode  =  12; 

int  pulsemode  =  1; 

int  countmode  =  2; 

float  delay; 
float  laser_freq; 
float  number_pulses; 
double  high_volt  =  3.50; 
double  low_volt  =  2.0; 
int  msec_pause; 
int  laser_half_count; 
char  retum_key[10]; 
static  char  *tname[]  =  { 

u  u 

» 

“Laser  Stop  Pulse”, 

“Laser  Pulse  Window”, 

“Laser  Shutter  Opening”, 

“Laser  Start  Gate”, 

“Camera  Shutter  Opening”, 

“Camera  Ready  Trigger”, 


u  u 


}; 

static  char  graph[]  =  (“trig. graph”}; 

double  start_gate_width  =  1000; 

double  stop_gate_width  =  1000; 

double  start_gate_delay,stop_gate_delay; 

double  pulse_window_width,  pulse_window_delay; 

double  camera_trigger_width  =  1400; 

double  wret.dret; 

int  camera_ready_pretrig; 

int  lshut_time_to_open  =  10000; 

int  high_pulse  =  1; 

int  retrig_rising_edge  =  3; 

mainO 

{ 

system(“  clear”); 
camera_ready_pretrig  =  10100; 
laser_freq  =  2000.; 

printfT  input  a  laser  frequency  :  “); 
scanfT%F,&laser_freq); 


while  Gaser  freq  >=  250  &&  laser  fireq  <=  10000) 

( 

mrclosallO; 

get_clocks_ready(); 

printfT  time  for  the  laser  shutter  to  open  %d  \n  “,lshut_time_to_open); 


/*  ALL  TIMES  ARE  IN  MICRO  SECONDS  */ 
start_gate_delay  =  camera_ready_pretrig  -  lshut_time_to_open; 
iflstart_gate_delay  <=  0) 

( 

printfT  Start  Gate  Delay  is  negative  !!  \n”); 
printfT  Increase  camera  pretrig  level  \n”); 
exit(l); 

} 

number_pulses  =  3; 


printfT  Delay  time  between  Camera  ready  pulse  and  Start  Pulse  :  %f  microsec 
\n”,start_gate_delay); 

printft"  Laser  External  Driver  Frequency  (square  wave)  :  %f  Hz  \n”,laser_freq); 

printfl"  Number  of  visible  pulses  requested  :  %f  \n”,number_pulses); 

printf(“  Time  it  takes  for  the  shutter  to  open  %d  \n  lshut_time_to_open); 


A-4  - 


/*  Make  sure  to  get  the  requred  number  of  pulses  */ 

/*  number  of  requested  pulses  plus  1  pre-ionization  pulse  and  */ 
/*  room  for  all  pulses  to  be  completed  */ 
number_pulses  =  number_pulses  +  1.05; 


set_clocks(); 

show_data(); 


printfT  CLOCKS  ARE  RUNNING  \n”); 


printfT  input  a  laser  frequency  :  “); 
scanft“%F,&laser_freq); 

/* 

printfT  Input  Camera  pretrig  :  “); 
scanfT%d”,&camera  ready_pretrig); 

*/ 

systemT  clear”); 

} 

/* 

printfT  Hit  Return  Key  to  Quit :  “); 
scanfT%c”,return_key); 

*/ 

mrclosalK); 

printfT  All  Devices  have  Been  Closed  \n”); 

} 


get_clocks  readyO 

{ 

printfT"  Get  the  clocks  assigned  and  ready  \n”); 
nclocks  =  4; 

for  ( i  =  0  ;  i  <  nclocks  ;  ++i) 

( 

clock_path[i]  =  -1; 
clk_da[i]  =  -1; 

) 

mropen(&clock_path[0],7dev/dacp0/clk3",0); 
mropen(&clock_path[l],7dev/dacp0/clk4",0); 
mropen(&clock_path[2]  ,7dev/dacp0/clk5",0); 
mropen(&clock_path[3],7dev/dacp0/clk6",0); 

mropen(&clk_da[0],  7dev/dacp0/clkda0”,  0); 
mropen(&clk_da[l],  “/dev/dacpO/clkdal”,  0); 
mropen(&clk_da[2],  7dev/dacp0/clkda2”,  0); 
mropen(&clk_da[3],  7dev/dacp0/clkda3",  0); 
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mrclkdis(nclocks,clock_path); 


printft"  set  the  trigger  level  to  %f  and  %f  volts  \n”, 

high_volt,  low_volt); 
mrclkdaset(clk_da[0],high_volt); 
mrclkdaset(clk_da[l],high_volt); 
mrclkdaset(clk_da[2],low_volt); 
mrclkdaset(clk_da[3],high_volt); 


1 

set_clocks() 

{ 

printft"  set  up  the  clocks  with  the  requested  values  \n”); 


/*  Start  Pulse  output  */ 

mroneshot(clock_path[0],0,start_gate_delay,&dret,0,start_gate_width, 

&wret,high_pulse,retrig_rising_edge); 

/*  Laser  Trigger  Gate  */ 

pulse_window_delay  =  start_gate_delay  +  lshut_time_to_open; 
pulse_window_width  =  number_pulses  /  laser_freq  *  1000000; 
mroneshot(clock_path[2],0,pulse_window_delay,&dret,0, 

pulse_window_width,&wret,high_pulse,retrig_rising_edge); 


/*  Put  out  the  laser  drive  frequency  ,  square  wave  */ 
laser_half_count  =  3000000.0  /  laser_freq 

mrclksetter(clock_path[3],sixmhz,laser_half_count,laser_half_count, 

highlevelgate,startlow,countinode,20); 


/*  Stop  Pulse  Output  */ 

stop_gate_delay  =  start_gate_delay  +  lshut_time_to_open  +  pulse_window_width 
mroneshot(clock_path[l],0,stop_gate_delay,&dret,0,stop_gate_width, 
&wret,high_pulse,retrig_rising_edge); 

printfT  Start  gate  delay  (usee) :  %f  \n”,start_gate_delay); 
printft"  Start  gate  width  (usee) :  %f  \n”,start_gate_width); 
printf(“  Pulse  Window  delay  (usee) :  %f  \n”,pulse_window_delay); 
printfl“  Pulse  Window  width  (usee) :  %f  \n”,pulse_window_width); 
printft"  Laser  half  count  :  %d  \n”,laser_half_count); 

printft"  Stop  gate  delay  (usee) :  %f  \n”,stop_gate_delay); 
printft"  Stop  gate  width  (usee) :  %f  \n”,stop_gate_width); 


mrclkarm(nclocks,clock_path); 


} 


show_data() 

{ 
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float  t[100]; 
float  v[  100]; 
float  start  =  0.0; 
float  offset ; 

float  cam_shut_open_time 
float  cam_shut_rise_time; 
float  restart; 

offset  =  100; 

restart  =  1000000.0/  25.0; 
cam_shut_rise_time  =  7000; 
cam_shut_open_time  =  1000; 


mpinit(gls); 

mpaxtype(gls,7,21); 

t[0]  =  -5000; 

/*  Camera  Trigger  */ 

v[0]  =  offset; 

t[l]  =  start; 

v[l]  =  offset; 

t[2]  =  t[l]; 

v[2]  =  offset  +  5; 

t[3]  =  start  +  camera_trigger_width; 

v[3]  =  offset  +  5; 

t[4]  =  t[3]; 

v[4]  =  offset; 

t[5]  =  start  +  restart; 

v[5]  =  offset; 

mplotsrcxf  gls,l,6,0,t,”F”,l,l,NTJLL,NULL); 
mplotsrcyf  gls,  l,6,0,v,”F”,  1, 1  ,NULL,NULL); 


offset  =  offset  -10; 

/*  Camera  Shutter  */ 
v[0]  =  offset; 

t[l]  =  start  +  camera_ready_pretrig  -  cam_shut_rise_time; 
v[l]  =  offset; 

t[2]  =  t[l]  +  cam_shut_rise_time; 
v[2]  =  offset  +  5; 

t[3]  =  t[2]  +  cam_shut_open_time; 
v[3]  =  offset  +  5; 

t[4]  =  t[3]  +  cam_shut_rise_time; 

v[4]  =  offset; 

t[5]  =  start  +  restart ; 

v[5]  =  offset; 

mplotsrcxf  gls, 2,6,0,t,”F”,  1,1, NULL, NULL); 
mplotsrcy(  gls,2,6,0,v,”F”,l,l,NULL,NULL); 


offset  =  offset  -10; 

/*  start  Gate  */ 
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v[0]  =  offset; 

t[l]  =  start  +  start_gate_delay; 

v[l]  =  offset; 

t[2]  =  t[l]; 

v[2]  =  offset  +  5; 

t[3]  =  t[2]  +  start_gate_width; 

v[3]  =  offset  +  5; 

t[4]  =  t[3]; 

v[4]  =  offset; 

t[5]  =  start  +  restart ; 

v[5]  =  offset; 

mplotsrcxC  gls,3,6,0,t,”F”,l,l,NULL,NULL); 
mplotsrcyC  gls,3,6,0,v,”F”,l,l,NULL,NULL); 

offset  =  offset  -10; 

/*  laser  Shutter  */ 
v[0]  =  offset; 

t[l]  =  start  +  start_gate_delay; 
v[l]  =  offset; 

t[2]  =  t[l]  +  lshut_time_to_open; 
v[2]  =  offset  +  5; 

t[3]  =  t[2]  +  pulse_window_width; 

v[3]  =  offset  +  5; 

t[4]  =  t[3]  +  lshut_time_to_open; 

v[4]  =  offset; 

t[5]  =  start  +  restart ; 

v[5]  =  offset; 

mplotsrcxC  gls,4,6,0,t,BF”,  1 , 1, NULL, NULL) ; 
mplotsrcyf  gls,4,6,0,v,”F”,  1,1,NULL,NULL); 

offset  =  offset  -10; 

/*  laser  Gate  */ 
v[0]  =  offset; 

t[l]  =  start  +  pulse_window_delay; 
v[l]  =  offset; 
t[2]  =  t[l]; 
v[2]  =  offset  +  5; 

t[3]  =  t[2]  +  pulse_window_width; 

v[3]  =  offset  +  5; 

t[4]  =  t[3]; 

v[4]  =  offset; 

t[5]  =  start  +  restart ; 

v[5]  =  offset; 

mplotsrcxC  gls,5 ,6,0,t,”FB,  1 , 1 , NULL,  NULL) ; 
mplotsrcyC  gls,5,6,0,v,”F”,l,l,NULL,NULL); 


offset  =  offset  -10; 

/*  stop  Gate  */ 
v[0]  =  offset; 

ttl]  =  start  +  stop_gate_delay; 
v[l]  =  offset; 
t[2]  =  t[l]; 
v[2]  =  offset  +  5; 
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t[3]  =  t[2]  +  stop_gate_width; 

v[3]  =  offset  +  5; 

t[4]  =  t[3]; 

v[4]  =  offset; 

t[5]  =  start  +  restart ; 

v[5]  =  offset; 

mplotsrcxf  gls,6,6,0,t,”F”,l,l,NULL,NULL); 
mplotsrcyf  gls,6,6,0,v,”F”,l,l,NULL,NULL); 

mpaxvalsfgls,  1, -5000.0, UNDEF, 5000.,- 1); 
mpaxvals(gls, 2, 40.0,1 10.0, 10. 0,-1); 

mptitle(gls,l,-l,-l,”  5  milliseconds  per  division”); 
mplabel(gls,2,7,l,0.0,-l,-l,tname); 


mpfilefgls,  graph,  1,2); 

mplot(gls,  0,1,0); 
mpend(gls); 


} 

transient_clock() 

{ 

/*  This  part  is  for  a  varing  laser  frequency  */ 

printfT  HIT  THE  RETURN  KEY  TO  START  RUN  ::  “); 

scanfT%c”,retum_key); 
printfT  countdown  starting  \n”); 

printfT  ...  5  ...  \n”); 
astpause(0, 1000); 
printfT"  ...  4  ...  \n”); 
astpause(0,1000); 
printfT"  ...  3  ...  \n”); 
astpause(0,1000); 
printfT"  ...  2  ...  \n”); 
astpause(0,1000); 
printfT"  ...  1 ...  \n”); 
astpause(0, 1000); 

printfT"  ...  going  with  first  clock...  \n”  ); 

msec  .pause  =  500; 
astpause(0,msec_pause); 


laser_freq  =  laser_freq  *2  ; 
set_clocks(); 

astpause(0,msec_pause); 


A-9/A-10 
Reverse  Blank 


Appendix  B 


FILM  PROCESSING  CONTROL  CODE:  RUNPIV.C 


B-l/B-2 
Reverse  Blank 


#include  <stdio.h> 

#include  <ctype.h> 

#include  <conio.h> 

#include  <math.h> 

int  i  j.k.l.m.n; 

int  start_frame,  number_frames; 
int  tid_number; 

char  command_string[40]; 

FILE  *fp,*fopen(); 

mainO 

( 

system(“cls”); 

printfl"  Welcome  to  the  Particle  Image  Velocimetry  Computing  System  \n”); 
printft"  The  system  will  compute  the  displacments  of  the  particles  \n\n”); 

printfl"  Make  sure  the  film  is  aligned  !  \n\n”); 

printfl"  Enter  the  Test  ID  number  :  “); 
scanfi“%d”,&tid_number); 

printfl"  Enter  the  First  frame  number  :  “); 
scanfi“%d”,&start_frame); 

printfi"  Enter  the  number  of  frames  :  "); 
scanfl“%d”,&number_frames); 


forfi  =  start_frame;  i  <  start_frame  +  number_frames;  i++) 

( 

system("cls”); 

printfl"  First  Digitize  the  picture  \n”); 
fp=fopen(“command.bat”,”w”); 

fprintfTfp,”  scanpic  ID%dFr%d.img  \n”,tid_number,i); 
fclose(fp); 

system(“command.bat”); 

printft"  Send  the  picture  to  cray  and  submit  the  job  \n”); 
fp=fopen(“command.bat”,”w”); 
fprintftfp,”  submit  ID%dFr%d  \n”,tid_number,i); 
fclose(fp); 

system(“command.bat”); 

printft"  Advance  the  picture  for  the  next  frame  \n”); 
system(“advance”); 


} 


) 


B-3/B-4 
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SCANNER  CONTROL  CODE:  SCANPIC.C 


C-l/C-2 
Reverse  Blank 


/*  Image  Scanning  program  by  Kenneth  LaPointe  Last  REV  2/21/92  */ 
/*  Reads  in  a  monochrome  image  from  Nikon  Film  Scanner  */ 

/*  compile  using  Microsoft  C,  "cl  scanpic.c  mcib.obj”  */ 


#include  <fcrtl.h> 

#include  <stdio.h> 

#include  <stdlib.h> 

#include  <string.h> 

#include  “decl.h” 

char  buffer[4600];  /*  buffer  array  to  transfer  data  using  GPIB  */ 
short  color[3200];  /*  storage  array  for  pixel  grey  scale  data  */ 
long  bufferjength;  /*  total  size  of  buffer  to  be  transfered  in  */ 
long  buffer_out;  /*  size  of  buffer  (bytes)  to  be  stored  in  file  */ 
int  header  =  4;  /*  size  of  header  in  GPIB  transfer  (words  )  */ 

int  etb  =  1;  /*  end  of  buffer  word  */ 

int  datajength;  /*  length  of  total  number  of  pixels  in  bytes  */ 

int  pitch;  /*  scanning  skip  ration  */ 

int  width;  /*  number  of  pixels  in  the  x  direction  to  scan  */ 

int  height;  /*  number  of  pixels  in  the  y  direction  to  scan  */ 

int  status;  /*  status  of  scanner  ,  busy/  not  busy  V 

int  ij,k,l,m,n,o;  /*  dummy  variables  */ 

char  file_name[30];  /*  array  to  hold  file  name  */ 

/*  FILE  VARIABLES  V 
int  datafiletn_write,n_read; 

main(argc.argv) 
int  argc; 
char  *argv[]; 

( 

/*  check  to  see  that  the  command  line  is  proper  */ 
iff  argc  ==111  argc  >=  3)( 

printff"  usage  :  scanpic  output_filename.img  \n”); 
exiffl); ) 

systemrcls”); 

datafile  =  open(argv[l],0_CREAT  I  O.TRUNC  1 0_BINARY  1 0_RDWR,0600); 
printff"  Data  file  :  %s  ;  has  been  opened  \n”, argv[l]); 

printff"  Make  the  GPIB  board  number  0  the  controller  in  charge  \n”); 
SendlFC(O); 

printff"  Clear  the  Scanner  ,  Device  address  is  1  \n”); 

DevClear(O.l); 

printff"  Set  the  scanner  for  monocrome  film  reading  a  negative  \n”); 

Send(0,  l,”MN\r”,3L,NLend); 

printff"  Auto  Focus  The  Film . \n”); 

Send(0,  l,"AF\r”,3L,NLend); 

WaitSRQ(0,  &status); 

printff"  Set  the  scanner  for  pixel  discharge  \n”); 
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Send(0,l,”SPOF\r”,5L,NLend); 

printft"  Do  a  prescan . \nB); 

Send(0,l,”PR\r”,3L,NLend); 

WaitSRQ(0,&status); 

printft"  Set  the  Exposure  time  Normal  =  40  :  \nB); 
Send(0,l,”AP40\rB,5L,NLend); 

printft"  Set  the  green  filter  for  black  and  white  \n”); 
Send(0,l,BFG\r”,3L,NLend); 

printft"  Requested  a  High  Res  Scan  \nB); 
high_res(); 


printft"  Return  scanner  to  zero  position  \nB); 

Send(0, 1  ,BMV0\  rB,4L,NLen  d); 

printft"  Close  the  Data  file  \nB); 
close(datafile); 

printft" . Finished  Scanning  Data . \nB); 

} 


high  rest) 

{ 

printft"  READ  HIGH  RES  4500  x  3000  \n  \nB); 
printft"  file  size  :  27  MB  \nB); 

/*  Start  at  250,1000  ,  End  at  4749,3999  V 
/*  area  of  box  size  4500  x  3000  */ 

/*  pitch  is  1  or  every  fifth  point  */ 

/*  number  of  data  points  900  x  600  V 
width  =  4500; 
height  =  3000; 
buffer_out  =  height  *  2; 

Send(0, 1,BBX  250, 1000, 1,1,4500, 3000, V\rB,28L,NLend); 
f*  I  */ 

r  1234567890123456789012345678901234567890  */ 

/*  buffer  length  is  size  of  header  plus  data  +  etb  */ 

/*  need  to  do  a  read  for  each  line  (pixel  array  )  of  data  */ 

/*  plus  one  last  blank  column  */ 
buffer_length  =  header  +  height  +  etb; 

printft"  READING  DATA  NOW . \nB); 

for(j  =  O  ;  j  <  widthj++) 

( 

/*  get  a  buffer  full  of  data  V 

Receive(0,l, buffer, buffer_length,STOPend); 

/*  Convert  the  character  data  to  integer  data  */ 

i  =  0; 

for(  k  =  header;  k  <  height  +  header;  k++) 
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( 

color[i]  =  bufferfk]; 

ifT  colorfi]  <  0)  color[i]  =  colorfi]  +  256; 
i++; 

) 

n_vvrite  =  write(datafile,  color,  buffer_out); 

) 

printfT"  Finished  with  data  \n”); 

/*  empty  the  last  blank  buffer  */ 

Receive(0,l, buffer, 5L.STO  Pend); 


I*  ======================================================  V 

low  res() 

{ 

printfl[“  READ  LOW  RESOLUTION  900  x  600  \n  \n"); 

printfT"  file  size  :  10  MB  \n”); 

/*  Start  at  250,1000  ,  End  at  4749,3999  V 
/*  area  of  box  size  4500  x  3000  */ 

/*  pitch  is  5  or  every  fifth  point  */ 

/*  number  of  data  points  900  x  600  */ 
width  =  900; 
height  =  600; 
buffer_out  =  height  *  2; 

Send(0,l,”BX  250,1000,5,5,900,600,V\r”,26L,NLend); 

/*  I  */ 

/*  1234567890123456789012345678901234567890  */ 

/*  buffer  length  is  size  of  header  plus  data  +  etb  */ 

/*  need  to  do  a  read  for  each  line  (pixel  array  )  of  data  */ 

/*  plus  one  last  blank  column  */ 
buffer_length  =  header  +  height  +  etb; 

printfT"  READING  DATA  NOW . \n”); 

for(j  =  0  ;  j  <  width  j++) 

{ 

/*  get  a  buffer  full  of  data  */ 

Receive(0, 1, buffer, buffer Jength.STOPend); 

/*  Convert  the  character  data  to  integer  data  V 

i  =  0; 

for(  k  =  header;  k  <  height  +  header;  k++) 

( 

colorfi]  =  buffer  [k]; 

ifT  colorfi]  <  0)  colorfi]  =  colorfi]  +  256; 
i++; 

) 

n_write  =  writefdata  file, color, buffer_out); 

) 

printfT"  Finished  with  data  \n”); 

/*  empty  the  last  blank  buffer  */ 

Receive(0, 1,  buffer,  5L.STO  Pend); 
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FILM  TRANSPORT  CONTROL  CODE:  ADVANCE.C 


D-l/D-2 
Reverse  Blank 


/*  Program  :  advance,  c  */ 

/*  Last  Rev  :  8/26/92  */ 

/*  Written  by  Ken  LaPointe  Code  8322,  NUWC  */ 
/*  */ 

/*  This  code  runs  the  film  movement  stages.  It  */ 

/’*'  will  move  the  film  one  frame  at  a  time  */ 

/*  */ 

/*  Compile  the  code  by  using  the  line  */ 

/*  cl  advance.c  pmc300.obj  */ 

/*  Make  sure  that  pmc300.obj  and  pmc300.h  are  */ 
/*  in  the  current  directory,  the  pmc  files  come  */ 

/*  from  the  Newport  software  */ 

#include  <stdio.h> 

#include  <ctype.h> 

#include  <conio.h> 

#include  <math.h> 

include  “PMC300.h” 

#define  x_axis  1 
#define  y_axis  2 
#define  T_axis  3 
#define  conv_mm  20000 
#define  conv_deg  1000 

float  frame_length  =  38.5; 

double  samfreq=  488.25; 
double  gain_x  =  50.0; 
double  pole_x  =  -0.1; 
double  zero_x  =  0.8; 
double  maxvel_x  =  0.40; 
double  accel_x  =  10.0; 
double  gain_y  =  50.0; 
double  pole_y  =  -0.1; 
double  zero_y  =  0.8; 
double  maxvel_y  =  0.40; 
double  accel_y  =  10.0; 
double  gain_r  =  10.0; 
double  pole_r  =  -0.1; 
double  zero_r  =  0.7; 
double  maxvel_r  =  8.00; 
double  accel_r  =  10.00; 

long  x_actj_act,r_act; 

int  i  j,k,l,m,n,o,p; 
int  frame_number; 
int  motion, stop; 

FILE  *fp,*fopen(); 

main() 
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( 

/***  starting  check  ********/ 
initialize_hardware(); 

/*******  move  film  ******/ 
move_film(); 

) 


1+*+*+*+++++:*+++++  subroutines  ***************/ 
y************************ **********************/ 


move_film() 

( 

/*  The  starting  position  is  close  to  the  NIKON  and  */ 

I*  away  from  the  film  */ 

/*  Now  Push  the  pins  into  the  film  */ 
x_act  =  5.0  *  conv_mm; 
set_com_pos(x_axis,x_act); 
check_motion(); 

/*  Move  the  Pins  to  Advance  the  film  one  frame  */ 
y_act  =  -  frame_length  *  conv_mm; 
set_com  _pos(y_axis,y_act); 
check_motion(); 

/*  Pull  the  pins  away  from  the  film  ,  the  zero  position  */ 
x_act  =  0; 

set_com_pos(x_axis,x_act); 

check_motion(); 

/*  Wind  the  Film  a  bit  in  degrees  */ 
r_act  =  -40  *  conv_deg; 
set_com  _pos(r_axis,r_act); 
check_motion(); 

/*  Return  the  pins  to  the  zero  position  for  the  next  advance  */ 
y_act  =  0; 

set_com_pos(y  axis,y_act); 

) 

check_motion() 

( 

motion  =  1; 
stop  =  0; 

while(motion  !=  stop) 

( 

movingO; 

} 

1 
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moving() 

{ 

long  counter I_old,cou.nter2_old,counter3_old; 
long  counterl_new,counter2_new,counter3_new; 
long  counter; 

int  x_motion,y_motion,r_motion; 


counterl_old  =  get_act_pos(x_axis); 
counter2_old  =  get_act_pos(y_axis); 
counter3_old  =  get_act_pos(r_axis); 

for(counter  =  0;  counter  <  100000;counter++) 
counter  =  counter; 

counter  l_new  =  get_act_pos(x_axis); 
counter2_new  =  get_act_pos(y_axis); 
counter3_new  =  get_act_pos(r_axis); 

x_motion  =  100; 
y_motion  =  10; 
r_motion  =  1; 

if(counterl_old  ==  counter  l_new) 
x_motion  =  0; 

if(counter2_old  ==  counter2_new) 
y_motion  =  0; 

ifTcounter3_old  ==  counter  3_new) 
r_motion  =  0; 

motion  =  x_motion  +  y_motion  +  r_motion; 

} 

initialize_hardware() 

{ 

mcJnitO; 

config_ports(l, 1,1,1); 

/***  set  up  x  axis  **+**/ 
reset(x_axis); 

set_sample_freq(x_axis,samfreq,0); 

set_zero(x_axis,zero_x); 

set_pole(x_axis,pole_x); 

set_gain(x_axis,gai  n_x); 

clr_act_pos(x_axi  s ) ; 

enter_ctl_mode(x_axis); 

set_prop_vel(x_axis,maxvel_x); 


/***  set  up  y  axis  ***++/ 
reset(y_axis); 
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set_sample_freq(y_axis,samfreq,0); 

set_zero(y_axis,zero  _y); 

set_pole(y_axis,pole  _y); 

set_gain(y_axis,gain_y); 

clr_act_pos(y_axis); 

enter_ctl_mode(y_axis); 

set_prop_vel(y_axis,maxvel_y); 

/***  set  up  rotary  axis  *****/ 
reset(r_axis); 

set_sample_freq(r_axis,samfreq,0); 

set_zero(r_axis,zero_r); 

set_pole{r_axis,pole_r); 

set_gain(r_axis,gain_r); 

clr_act_pos(r_axis); 

enter_ctl_mode(r_axis); 

set_prop_vel(r  axis.maxvel  r); 

} 
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CRAY  JOB  CONTROL  CODE:  SUBMIT.BAT 


E-l/E-2 
Reverse  Blank 


©echo  off 

echo  Script  to  send  the  Image  file  to  the  CRAY 
echo  and  compute  the  displacement  vectors 
echo  the  data  will  be  sent  to  the  PI83a  computer 


echo  Make  the  Job  script 
echo  rm  %l.e*  %l.o* 
echo  cd  /usr/tmp/PIV 
echo  rm  %l.act 
echo  rm  %l.dsp 
echo  ja 

echo  fftpiv  %l.img  %l.dsp 

echoja  -st 

echo  %l.snd 

echo  rm  %l.img 

echo  rm  %l.job 

echo  rm  %l.snd 


>  %l.job 
»  %l.job 
»  %l.job 
»  %l.job 
»  %l.job 

»  %l.job 
»  %l.job 
»  %l.job 
»  %l.job 
»  %l.job 
»  %l.job 


echo  Make  The  Cray  Send  the  info  back  to  the  SGI 
copy  send.h  %l.snd 
echo  put  %l.dsp  »  %l.snd 

echo  put  %l.act  »  %l.snd 

echo  cd  »  %l.snd 

echo  mput  %1*  »  %l.snd 

echo  EOF  »  %l.snd 


echo  SET  UP  THE  PC  TO  SEND  THE  INFO  TO  THE  CRAY 
copy  C:\ENET\cray.log  %l.put 


echo  prompt 
echo  verbose 
echo  ascii 

echo  cd  /usr/tmp/PIV 
echo  put  %l.job 
echo  put%l.snd 
echo  binary 
echo  put  %l.img 
echo  bye 


»  %l.put 
»  %l.put 
»  %l.put 

»  %l.put 
»  %l.put 
»  %l.put 
»  %l.put 
»  %l.put 
»  %l.put 


C:\PATHWAY\ftp  <  %l.put 


echo  Make  the  %l.snd  file  exacutable 

rsh  Cray  chmod  777  /usr/tmp/PIV/%l.snd 


echo . 
echo . 


echo  ********  check  on  files  sent  ******* 

C:\PATHWAY\rsh  cray  Is  -al  /usr/tmp/PIV 

echo  Submit  the  batch  Job  to  the  Cray 

rsh  cray  qsub  -IT  00:08:00  -lm  2.5MW  -r  %1  /usr/tm p/PI V/%  1  .j ob 
rsh  cray  qstat  -ba 


echo  remove  the  files  on  the  PC 
del  %l.img 
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del  %l.snd 
del  %l.put 
del  %l.job 
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AUTOCORRELATION  PIV  CODE:  AUTOCORE.C 


F-l/F-2 
Reverse  Blank 


/**  PARTICLE  IMAGE  VELOCEMETRY  PROGRAM  **/ 
/**  using  the  autorcorrelation  method  **/ 

/**  By  :  Kenneth  LaPointe  **/ 

/**  NUSC,  Code  8322  **/ 

/**  Last  rev  as  of  AUG  17  1992  **/ 

/**  This  program  reads  an  image  **/ 

/**  files  to  use  an  auto  correlation  program  **/ 

/**  to  compute  the  displacemnet  between  **/ 

/**  particles  **/ 

/**  compile  on  sgi  as  :  cc  -03  autocore.c  -Dsgi  -lm  -o  autocore  */ 
/**  compile  on  cray  as  :  cc  -03  autocore.c  -Dcray  -lm  -o  autocore 

/**  Input  for  the  program  is  *  V 

/**  autocore  filename.img  filename. dsp  filename. up 

/**  the  .img  file  is  the  input  file  scanned  by  the  NIKON  scanner 
/**  the  .dsp  file  is  the  output  file  computed  by  the  program 
/**  the  .up  is  the  setup  file  for  the  code  */ 

#include  <stdio.h> 

#include  <math.h> 

#include  <unistd.h> 

#include  <sys/types.h> 

int  true  =  1; 
int  false  =  0; 
int  i  j,k,l,m,n,o,p; 
float  a,b,c,d,e,f,g; 

int  junk; 
int  nbuff; 
int  buffer_bytes; 
int  buf_loop,x_offset; 
int  pic_width; 
int  win_rel; 


*/ 


V 

*/ 

V 


short  rgbvector[3]; 
long  datastore; 


intxwl,ywl;  /*  lower  left  corner  of  2  D  array  is  xwl,ywl  */ 

int  xw2,yw2;  /*  upper  right  comer  of  2  D  array  is  xw2,yw2  */ 

int  csx,csy;  /*  copy  shift  value  in  x  and  y,  is  the  amount  the  copy  is  moved  */ 

int  start_csx,start_csy;  /*  amount  of  initial  movment  of  the  copy  */ 

int  end_csx,end_csy;  /*  move  the  copy  to  this  point  */ 

int  hold_csy;  /*  hold  buffer  varible  for  the  csy  value  */ 

inthold_csx; 

int  isx.isy;  /*  2D  image  array  index  start  point  in  x  and  y  */ 

int  iey,iex;  /*  2D  iamge  array  index  end  point  in  x  and  y  */ 


int  iix.iiy;  /*  2D  image  array  index  x,y  V 

int  cix,ciy;  /*  2D  copy  array  index  x,y,  */ 

int  range_x,range_y; 

int  xw_width,yw_width; 
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/*  for  the  correlatons,  there  are  three  positions,  */ 

int  pz  =  2;  /*  number  2  is  the  zero  position  */ 

int  pn  =  1;  f*  number  1  is  the  minus  position  */ 

intpp  =  0;  /*  number  0  is  the  positive  position  */ 

int  xp  =  0;  /*  and  the  x  and  y  values  */ 

intyp  =  1; 

int  disp_peak[3][2]; 

int  scan_area; 

float  corr_hist[2][200][200]; 

float  corr_peak[3]; 


int  n_windows,win; 
int  max_frames,run_number; 
int  x_win_step,y_win_5tep; 
float  total; 

int  max_x_disp,max_y_disp; 

char  file_name[20];  /*  file  name  character  storage  */ 

char  last_data_file[20]; 


int  image_file;  /*  file  pointer  for  image  file  */ 

int  height, width; 

int  height; 

char  old_disp_name[50]; 
char  new_disp_name[50]; 
int  shade; 

int  n_wind_sec_x,  n_wind_sec_y; 

int  frame_number; 

int  dia_part; 

int  xwin,ywin; 

float  skipx,  skipy; 

int  range_x_guess,  range_y_guess; 

int  cs_csx,ce_csx,cs_csy,ce_csy; 
int  draw_data; 
int  shift; 


float  x_sum,y_sum,x_weight,y_weight; 
float  x_sum_total,y_sum_total; 

int  shade_limit; 

#ifdef  cray 

unsigned  short  image_data[512][3200J; 
long  buffer[1024]; 
unsigned  long  sl,s2,s3,s4; 
off_t  seek_error,image_offset_bytes; 
int  max_array_width  =  500; 
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int  max_array_height  =  3000; 
int  peak_value[500][2); 
float  avg_value[500][2]; 
float  zero_corr_ratio[500]; 

#endif 

#ifdef  sgi 

unsigned  short  image_data[512][3200]; 

short  buffer[3 100]; 

off_t  seek_error,image_offset_bytes; 

int  max_array_width  =  500; 

int  max_array_h  eight  =  3000; 

int  peak_value[500][2]; 

float  avg_value[500][2]; 

float  zero_corr_ratio[500]; 

#endif 

int  fd,n_read,n_write, bytes;  /*  file  descriptors  */ 

FILE  *file_dat,*look,*sfile,1''new_file,*old_file,*fopen(); 

j*  I*************************************************************** 


main(argc,argv) 
int  argc; 
char  *argv[]; 

{ 

/*  =====================  REQUIRED  INPUTS  ==================  */ 

/*  A  35  mm  negative  is  36mm  x  26  mm  ,  with  the  Nikon  scanner  at  a  pitch  of  1  */ 

/*  the  negative  is  digitized  to  125  pixels/mm  or  4500  x  3000  pixels  */ 

/*  For  use  on  the  Cray  the  size  of  the  buffer  and  picture  must  be  an  integer  */ 
/*  multiple  of  4  */ 

/**  read  in  the  set  up  file  */ 
sfile  =  fopen(argv[3],V); 

fscanftsfile,”%d”,&height); 

fscanflsfile,”%d”,&width); 

fscanfisfile,”%d”,&xw_width); 

fscan  ft  sfile, ”%d”,&yw_width); 

fscanflsfile,”%d”,&shade_limit); 

fscanfisfileI”%d”,&max_x_disp); 

fscanflsfile,"%d”,&max_y_disp); 

fscanfisfile,”%d”,&range_x  _guess); 

fscanfisfile,”%d”,&range_y  _guess); 

fscanfisfile1”%d”,&dia_part); 

fclose(sfile); 

/* 

Standard  settings  are 
height  =  3000; 
width  -  4500; 
xw_width  =  80; 


F-5  . 


yw_width  =  80; 
shadejimit  =  10; 
max_x_disp  =  40; 
max_y_disp  =  40; 
range_x  _guess  =  30; 
range_y  _guess  =  30; 
dia_part  =  12; 

*/ 


printfT  height  %d  \n", height); 
printfT  width  %d  \n”, width); 
printfT  xw  width  %d  \n”,xw_width); 
printfT  yw  width  %d  \n”,yw_width); 
printfT  shade  limit  %d  \n”, shadejimit); 
printfT  max  x  disp  %d  \n”,max_x_disp); 
printft"  max  y  disp  %d  \n”,max_y_disp); 
printfT  range  x  %d  \n",range_x_guess); 
printfT“  range  y  %d  \n",range_y_guess); 
printfT  dia  part  %d  \n”,dia_part); 


iflheight  >  max_array_height ) 

l 

printfT"  height  too  big  \n”); 
exit(l); 

) 


#ifdef  sgi 

ifT width  >  max  array_width  ) 

( 

pic_width  =  width; 

width  =  max_array_width  -  max  x  disp  -  range_x_guess; 

) 

#endif 
tifdef  cray 

iff  width  >  max_array_width) 

{ 

pic_width  =  width; 
width  =  xw_  width; 

) 


#endif 


/*  for  the  each  input  buffer  */ 
n_wind_sec _y  =  height /yw_width; 

n_wind_sec_x  =  width  /  xw_width;  I*  number  of  windows  in  the  x  direction  */ 
n_windows  =  n_wind_sec_x  *  n_wind_sec_y^*  total  number  of  windows  and  vectors  *1 
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nbuff  =  (pic_width  /  xw_width)/  n_wind_sec_x; 

printff“  max  disp  %d  \n”,max_x_disp); 

printfT  range  x  guess  %d  \n”,range_x_guess); 
printft"  dia  part  %d  \n”,dia_part); 

printff“  n  sec  y  %d  \n”,n_wind_sec_y); 
printfT  n  sec  x  %d  \n”,n_wind_sec_x); 
printfl"  n  windows  %d  \n”,n_windows); 

printfT"  number  of  buffers  %d  \n”, nbuff); 


/***  open  the  Image  data  file  ***/ 

image_file  =  open(argv[l],0); 

/***  open  the  displacment  data  file  **/ 
new_file  =  fopen(argv[2],”w”); 
fprintflnew_file,”%d  \n",run_number); 
fprintfTnew_file,”%d  \ n”,frame_n umber); 
fprintflnew_file,”%d  \n”, width  *  nbuff); 
fprintflnew_file,”%d  \n”, height); 
fprintfinew_file,”%d  \n”,n_wind_sec_x  *  nbuff); 
fprintf(new_file,”%d  \n”,n_wind_sec_y); 


x_offset  =  0; 

#ifdef  cray 

/*  compute  the  number  of  buffers  */ 
buffer_bytes  =  height  *  2; 

for  (bufjoop  =  0;  buf  loop  <  nbuff;  bufJoop++) 

{ 

/*+**+***  for  testing  *******/ 

image_offset_bytes  =  buf_loop  *  xw_width  *  n_wind_sec_x  *  buffer_bytes; 

get_image_file(); 

pivO; 

dump_disp(); 


#endif 
#ifdef  sgi 

/*  compute  the  number  of  buffers  */ 
bufferjbytes  -  height  *  2; 

for  (bufjoop  =  0;  bufjoop  <  nbuff;  buf_loop++) 

( 

image_offset_bytes  =  bufjoop  *  xw_width  *  n_wind_sec_x  *  buffer_bytes; 

get_image_file(); 

pivO; 

dump_disp(); 

) 
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#endif 


/*  close  the  files  *1 
fclose(new_file); 
closed  mage  Jile); 

) 


#ifdef  cray 
get_image  file() 

( 

int  endjoop; 

/*  Read  in  all  of  the  data  line  by  line  */ 

/*  the  data  is  read  into  the  array  "  buffer”  */ 

/*  Since  the  data  went  from  a  16  bit  machine  (PC)  to  a  64  bit  machine  (CRAY)  */ 
/*  four  data  points  are  in  one  CRAY  word  */ 

/* 

printfi"  buf  loop  %d  width  %d  buf  bytes  %d  image  offset  %ld  xoffset  %d  \n”, 
bufjoop, width, bufferJ>ytes,image_offsetJ)ytes,x_offset); 

*/ 


seek_error  =  lseek(image_file,image_offset_bytes,SEEK_SET); 

end_loop  =  2  *  width; 

ifi  bufjoop  +  1  ==  nbuff) 

end_loop  =  xw_width  *  n_wind_sec_x; 

for  ( j  =  0  ;  j  <  end  loop  ;  ++j) 

{ 

m  =  0; 

n_read  =  read(image_file .buffer, height  *  2); 
for(  k  =  0  ;  k  <  height/4;  k++) 

{ 

image_data(j][m]  =  image_data[j](m+l]  =  image_data|J][m+2]  =  image_data[j][m+3]  =  0; 
si  =  buffer[k]  ; 
image_data(j][m  ]  =  si  »  56; 

if(image_data(j][m]  <  shadejimit  I  I  image_data[j][m]  >  255)  image_data(j)[m]  =  0; 

s2  =  buffer[k]  «  16; 
image_data[j][m+l]  =  s2  »  56; 

iflimage_data[j][m+l]  <  shade_limit  I  I  image_data(j][m+l]  >  255)  image_data[j][m+l]  =  0; 

s3  =  bufferfk]  «  32; 
image_data!j][m+2]  =  s3  »  56; 

ifiimage_data[j][m+2]  <  shadejimit  I  I  image_data(j][m+2]  >  255)  image_data[jKni+2]  =  0; 

s4  =  buffer[k]  «  48; 
image_data[j][m+3]  =  s4  »  56; 


F-8  • 


iffimage_data[j][in+3]  <  shadejimit  I  I  image_data[j][m+3]  >  255)  image_data[j][m+3]  =  0; 


m  +=4; 

) 


) 


) 

#endif 


#ifdef  sgi 
getjmage  file() 

( 

intendjoop; 

printff“  reading  data  An"); 

/*  Read  in  all  of  the  data  line  by  line  */ 

/*  the  data  is  Tead  into  the  array  “  buffer”  */ 

/*  Since  the  data  went  from  a  16  bit  machine  (PC)  to  a  32  bit  machine  (SGI)  */ 

/*  each  element  of  the  array  is  bit  shifted  by  8  bits  */ 

/*  The  data  on  the  PC  is  l’s  complement,  the  data  on  the  SGI  is  2’s  */ 

/*  complement,  so  if  the  data  is  less  than  0  add  256  V 


seek_error  =  lseek(image_file,image_offset_bytes,SEEK_SET); 

iffseek_error  ==  -1) 
exit(l); 

endjoop  =  width; 

iff  bufjoop  -1  ==  nbuff) 

endjoop  =  xw_width  *  n_wind_sec_x; 
for  ( j  =  0  ;  j  <  end  loop  ;  ++j) 

( 

n_read  =  read(image_file, buffer, height  *  2); 
for(  k  =  0  ;  k  <  height;  ++k) 

( 

shade  =  buffer[k]  »  8; 
iff  shade  <  0  )  shade  +=  256; 
iff  shade  <  shadejimit  )  shade  =  0; 
image_data[j][k]  =  shade; 

) 


} 

#endif 


dump_disp() 
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{ 

win  =  1; 

for  (x_win_step  =  0  ;  x_win_step  <  n_wind_sec  x;  x_win_step++) 

{ 

xwl  =  x_win_step  *  xw_width  +  x_offset; 
xw2  =  xwl  +  xw_width; 

for  (y_win_step  =  0  ;  y_win_step  <  n  wind_sec_y;  y_win_step++) 

( 

ywl  =  y_win_step  *  yw_width; 
yw2  =  ywl  +  yw_width; 

fprintfTnew_file,,’%d  %d  %d  %d  %d  %d  %d  %f  %f  %f\n”,win_rel,xwl,ywl,xw2,yw2, 
peak_value[win][xp],peak_value[win][yp],avg_value[win][xp], 
avg_value[win][yp], 
zero_corr_ratio[win]); 

printfl“%d  %d  %d  %d  %d  %d  %d  %f  %f  \n’,,win_rel,xwl,ywl,xw2,yw2, 
peak_value[win](xp],peak_value[win][yp],avg_value[win][xp], 
avg_value[win][yp]); 

win++; 

win_rel++; 

) 

) 

x_offset  =  xw2; 

) 


piv() 

( 

xwl  =  0;  /*  starting  point  of  window  in  x  direction  (in  pixels)  V 

ywl  =  0;  /*  y  */ 

win  =  1;  f*  first  window  number  */ 

/*  in  This  routine  loop  through  each  window,  get  the  two  comers  of  the  window  */ 

/*  use  the  maximumm  estamated  displacment  to  compute  the  auto-correlation  */ 

for  (x_win_step  =  0  ;  x  win_step  <  n_wind  sec_x;  x_win_step++) 

{ 

xwl  =  x_win_step  *  xw_width; 
xw2  =  xwl  +  xw_width; 
range_x  =  max_x_disp; 
range_y  =  max_y_disp; 

for  (y_win_step  =  0  ;  y_win_step  <  n  wind_sec _y;  y  win  step++) 

{ 

ywl  =  y_win_step  *  yw_width; 
yw2  =  ywl  +  yw_width; 

printfl“  ******  win  %d  xwl  %d  xw2  %d  ywl  %d  yw2  %d  ***  \n", 
win,xwl,xw2,ywlyw2); 

auto_corr(); 
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printft“%d  %d  %d  %d  %d  %d  %d  %f  %f  \n",win_rel,xwl,ywl,xw2,yw2, 
peak_value[win]Lxpl,peak_value[win][yp],avg_value[win][xp], 
avg_valuefwin][yp]); 

range_x  -  range_x_guess; 
range_y  =  range_y _guess; 
win++; 

) 


) 


auto  corr() 

{ 

/*  get  the  best  possible  autocorrelation  value,  the  zero  position  */ 
perform_zero_corr(); 

/*  Scan  through  the  positive  x  and  y  dispacments  */ 
perform_pos_corr(); 

/*  Scan  through  the  positive  x  and  negative  y  dispacments  */ 
perform_neg_corr(); 


printfT  POS  DISPLACMENT  DATA  max  value  %.2f  at :  x  %d  y  %d  \n”, 
corr_peak[pp],disp_peak[pp][xp],disp_peak[pp][yp]); 


printft“  NEG  DISPLACMENT  DATA  max  value  %.2f  at :  x  %d  y  %d  \n", 
corr_peak[pn],disp_peak[pn][xp],disp_peak[pn][yp]); 


avg_value[win][xp]  =  0; 
avg_value[win][yp]  =  0; 

peak_value[win][xp]  =  o; 

peak_valu.e[win][yp]  =  0; 


ift  corr_peaktpp]  >  corr_peak[pn] ) 

{ 

/*  displacment  is  positive  y*/ 
scan_area  =  pp; 

zero_corr_ratio[win]  =  corr_peak[pp]  /corr_peak[pz]; 
scan  histO; 

) 

else 

( 

/*  displacment  is  negative  y*/ 
scan_area  =  pn; 

zero_corr_ratio[win]  =  corr_peak[pn]  /corr_peak[pz]; 
scan_hist(); 
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avg_value[win](yp]  =  -  avg_value[win][yp]; 

) 


perform  zero_corr() 

{ 

/*  upper  right  quad  */ 
corr_peak[pz]  =  0.0001; 
disp_peak[pz][xp]  =  0; 
disp_peak[pz]typ]  =  0; 

/**  Copy  is  placed  on  top  of  the  image  ,  zero  position**/ 

total  =  0; 

isy  =  y  w  1  ; 

iey  =  yw2  ; 

isx  =  xwl  ; 

iex  =  xw2  ; 

/**  perform  auto  correleation  for  this  copy  placecent  */ 
ciy  =  isy; 

for(  iiy  =  isy;  iiy  <  iey;  iiy++) 

( 

cix  =  isx; 

for(  iix  =  isx;  iix  <  iex;  iix++) 

l 

total  +=  image_data[iix][iiy]  *  image_data[cix][ciy]; 
cix++; 

) 

ciy++; 

} 

corr_peak[pz]  =  total; 

) 

perform_pos  corr() 

( 

/*  upper  right  quad  +x  ,  +y  movment  of  the  copy*/ 
corr_peak[pp]  =  0; 
disp_peak[pp][xp]  =  0; 
disp_peak[pp][yp]  =  0; 

printfT  Set  up  the  positive  corr  parameters  \n”); 

/*  Set  up  the  movment  positions  of  the  copy  on  the  original  */ 

/*  make  sure  the  copy  does  not  do  the  zero  position  and  */ 

/*  does  not  go  outside  the  bounds  of  the  array  memory  locations  */ 
/*  set  the  dispacment  starting  and  stopping  points  based  on  a  best  */ 
/*  guess  and  some  resonable  range  estimate  */ 

start_csy  =  abs(peak_value[win  -l][yp])  -  range_y  ; 
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ifTstart_csy  <  0  )  start_csy  =  0; 


end_csy  =  abs(peak_value[win  -l][yp])  +  range_y; 
iflend_csy  >  height  )end_csy  =  height; 

start_csx  =  abs(peak_value[win  -l][xp])  -  range_x; 
iftstart_csx  <  0  )  start_csx  =  0; 

end_csx  =  abs(peak_value[win  -l][xp])  +  range_x; 
iflend_csx  >  width)  end_csx  =  width; 

/*  Remember  not  to  do  the  zero  position,  The  zero  position  is  actually  */ 

/*  the  same  size  as  a  particle  is  in  pixels  */ 

/*  If  the  copy  is  moved  over  the  0,0  position  */ 

/ *  do  the  correlation  in  two  steps  */ 

/*  to  the  left  and  top  of  the  partilcle  and  then  all  area  above  it  */ 

printfT  requested  positive  corr  start  csy  %d  end  %d  start  csx  %d  end  %d  \n”, 
start_csy,end_csy,start_csx,end_csx); 

if[start_csy  <-  dia_part  &&  start_csx  <=  dia_part ) 

{ 

hold_csy  =  start_csy; 
start_csy  =  dia_part; 
hold_csx  =  end_csx; 
end_csx  =  dia_part; 
pos_copy_shift() ; 

start_csy  =  hold_csy; 
start_csx  =  dia_part; 
end_csx  =  hold_csx; 
pos_copy_shift(); 

} 

else 

( 

pos_copy_shift(); 

} 

) 


pos  copy_shift() 

{ 

printft"  actual  positive  corr  start  csy  %d  end  %d  start  csx  %d  end  %d  \n”, 
start_csy,end_csy,start_csx,end_csx); 

/**  shift  copy  in  x  and  y  **/ 

/*  Place  the  copy  in  the  requested  starting  position  */ 

for(csy  =  start_csy;  csy  <  end_csy;  csy++) 
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{ 

isy  =  ywl  +  csy  ; 
iey  =  yw2  +  csy  ; 

for(csx  =  start_csx;  csx  <  end_csx;  csx++) 

( 

total  =  0; 
isx  =  xwl  +  csx  ; 
iex  =  xw2  +  csx  ; 
ciy  =  ywl; 

for(  iiy  =  isy;  iiy  <  iey;  iiy++) 

( 

cix  =  xwl; 

for(  iix  =  isx;  iix  <  iex;  iix++) 

{ 


total  +=  image_data[iix][iiy]  *  image_data[cix][ciy]; 

cix++; 

} 

ciy++; 

) 


/* 

if(win_rel  ==  60  &&  win  ==  41) 

printf(“  %d  %d  total  value  is  %f  \n”,csx,csy, total); 

*/ 


corr_hist[pp][csx][csy]  =  total; 
if( total  >  corr_peak[0j) 

{ 

corr_peak[pp]  =  total; 
disp_peak[pp][xp]  =  csx; 
disp_peak[pp][yp]  =  csy; 
) 

)  /*  end  loop  csx  V 
)  /*  end  loop  csy  */ 

printfT  Finished  pos  auto  corr  \nB); 

) 


perform  neg_corr() 

( 

/*  lower  right  quad  */ 

corr_peak[pn]  =  0; 

disp_peak[pn][xp]  =  0; 

disp_peak[pn][yp]  =  0; 

printfT  Set  up  the  neg  corr  parameters  \n"); 

start_csy  =  -abs(peak_value[win  -l][yp])  -  range_y  ; 
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iflstart_csy  <  -  height )  start_csy  =  -  height; 

end_csy  =  -abs(peak_value[win  -l][yp])  +  range_y; 
iftend_csy  >=  0  )end_csy  =  0; 

start_csx  =  abs(peak_value[win  -l][xp])  -  range_x; 
iftstart_csx  <  0  )  start_csx  =  0; 

end_csx  =  abs(peak_value[win  -l][xp])  +  range_x; 
ifl[end_csx  >  width)  end_csx  =  width; 

printf(“  requested  negative  corr  start  csy  %d  end  %d  start  csx  %d  end  %d  \n”, 
start_csy,end_csy,start_csx,end_csx); 

iflend  csy  >=  -dia_part  &&  start_csx  <=  dia_part ) 

{ 

hold_csy  =  end_csy; 
end_csy  =  -  dia_part; 
hold_csx  =  end_csx; 
end_csx  =  dia_part; 
neg_copy_shift(); 

end_csy  =  hold_csy; 
start_csx  =  dia_part; 
end_csx  =  hold_csx; 

neg_copy_shift(); 

) 

else 

( 

neg_copy_shift(); 

} 


neg_copy_shift() 

{ 

printfT"  actual  negative  corr  start  csy  %d  end  %d  start  csx  %d  end  %d  \n”, 
start_csy,end_csy,start_csx,end_csx); 

/**  shift  copy  in  x  and  y  **/ 

/*  rember  to  use  the  absolut  values  for  the  array  index  */ 
foKcsy  =  start_csy;  csy  <  end_csy;  csy++) 

( 

isy  =  ywl  ; 

iey  =  yw2  -  abs(  csy ); 

forfcsx  =  start_csx;  csx  <  end_csx;  csx++) 

( 

total  =  0; 

isx  =  xwl  +  csx  ; 

iex  =  xw2  +  csx  ; 

ciy  =  ywl  +  abs(csy); 

for(  iiy  =  isy;  iiy  <  iey;  iiy++) 
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{ 

cix  =  xwl; 

for(  iix  =  isx;  iix  <  iex;  iix++) 

{ 


total  +=  image_data[iix][iiy]  *  image_data[cix][ciy]; 
cix++; 

} 

ciy++; 

) 

corr_hist[pn][csx][abs(csy)]  =  total; 
ifl  total  >  corr  peaktpn]) 

( 

corr_peak[pn]  =  total; 
disp_peak[pn][xp]  =  csx; 
disp_peak[pn][yp]  =  csy; 


) 

printfl“  Finished  neg  auto  corr  \n”); 

) 


sean_hist() 

( 

/*  The  histogram  yields  a  peak  which  can  be  used  to  compute  the  displacment  vector  */ 
/*  A  better  vector  would  be  obtained  by  computing  the  average  displacement  around  *1 
/*  the  peak.  Compute  the  average  by  getting  the  weighted  value  in  each  dimension  */ 

/*  Scan  the  correlation  map  and  compute  the  average  value  around  the  peak  */ 

x_sum  =  0; 
y_sum  =  0; 
x_sum_total  =  1; 
y_sum_total  =  1; 
x_weight  =  0; 
y_weight  =  0; 

cs_csx  =  disp_peak[scan_area][xp]  -  4; 
if  (cs_csx  <  0)  cs_csx  =  0; 

ce_csx  =  disp_peak[scan_area][xp]  +  4; 

cs_csy  =  disp_peak[scan_area][yp]  -  4; 
ce_csy  =  disp_peak[scan_area][yp]  +  4; 

/*  compress  the  x  axis  to  one  dimension  */ 

/*  and  weight  the  displacment  by  the  vauel  of  the  peak  */ 

/*  compute  the  average  displacment  in  the  x  direction  V 

for(csx  =  cs  csx;  csx  <=  ce  csx;  csx++) 

{ 

for(csy  =  cs_csy;  csy  <=  ce_csy;  csy++) 
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{ 

x  sum  +=  corr_hist[scan_area][csx][abs(csy)]; 

) 

x_weight  +=  csx  *  x_sum; 
x_sum_total  +=  x_sum; 
x_sum  =  0; 

} 

avg_value[win][xp]  =  x_weight/x_sum_total; 


/*  compress  the  y  axis  to  one  dimension  */ 

for<csy  =  cs_csy;  csy  <=  ce_csy;  csy++) 

( 

for(csx  =  cs_csx;  csx  <=  ce_csx;  csx++) 

( 

y  sum  +=  corr_hist[scan_area][csx][abs(csy)]; 

) 

y_weight  +=  csy  *  y_sum; 
y_sum_total  +=  y_sum; 
y_sum  =  0; 

} 

avg_value[win][yp]  =  y_weight/y_sum_total; 

peak_value[win][xp]  =  disp_peak[scan_area][xp]; 
peak_value[win][yp]  =  disp_peak[scan_area][yp]; 

} 
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Reverse  Blank 


Appendix  G 


2D  FFT  PIV  CODE  FOR  THE  SGI:  SGI_FFTPIV.C 


G-l/G-2 
Reverse  Blank 


/**  PARTICLE  IMAGE  VELOCEMETRY  PROGRAM  **/ 

/**  By :  Kenneth  LaPointe  **/ 

/**  NUSC,  Code  8322  **/ 

/**  Last  rev  as  of  May  5  1992  **/ 

/**  This  program  reads  in  data  files  and  image  **/ 

/**  files  to  use  an  2D  FFT  program  **/ 

/**  to  compute  the  displacemnet  between  **/ 

/**  particles  **/ 

/**  compile  on  sgi  as  :  cc  -03  sgijftpiv.c  -lm  -o  sgi_fFtpiv  */ 

/*  To  see  the  IMAGE,  and  FFTs  on  the  screen  compile  with  */ 

/**  compile  on  sgi  as  :  cc  -03  sgi_fftpiv.c  -Dshow  -lgl  -lm  -o  sgi_fftpiv  */ 

#include  <stdio.h> 

#include  <math.h> 

#include  <unistd.h> 

#include  <sys/types.h> 

#ifdef  show 
#in elude  <gl/gl.h> 

#endif 

int  try; 

int  end_buflfer; 

int  fftjoop  ; 

int  fft_pass; 

int  true  =  1; 

int  false  =  0; 

intij.k.l.m.n.o.p; 

int  ajl,aj2,aj3,alcl,ak21ak3; 

float  a,b,c,d,e,f,g; 

float  cp; 

double  square,  magnitude; 

float  fft_peak; 
float  f!t_ratio; 
float  cen_ratio; 
float  dc_com; 

float  first_peak_value; 
int  first_peak_x; 
int  first_peak_y; 

float  second_peak_value; 
int  second_peak_x; 
int  second_peak_y; 
int  aindex; 

int  junk; 

int  nbuff; 

int  buffer_bytes; 

int  buf_loop,x_offset,start_buffer; 
int  pic_width ; 
int  win_rel; 

short  rgbvector[3]; 
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long  datastore; 


int  xwl,ywl;  /*  lower  left  corner  of  2  D  array  is  xwl.ywl  */ 

int  xw2,yw2;  /*  upper  right  comer  of  2  D  array  is  xw2,yw2  */ 

int  xw_width,yw_width; 

int  n_windows,win; 
int  max_frames,run_number; 
int  x_win_step,y_win_step; 
int  total; 

int  xwlout,xw2out; 


int  max_x_disp,max  _y_disp; 

char  file_name[20];  /*  file  name  character  storage  */ 

char  last_data_file[20]; 


int  image_file; 

int  image.height, image,  width; 

char  old_disp_name[50]; 
char  new_disp_name[50]; 
int  shade; 

int  n_wind_sec_x,  n_wind_sec_y; 

int  frame.n  umber; 

int  dia_part; 

int  xwin,ywin; 

float  skipx,  skipy; 

float  cp; 
int  draw.data; 
int  shift; 

int  peak_value_x; 
int  peak_value_y; 
float  avg_value_x; 
float  avg_  value _y; 

float  x_sum,y_sum,x_weight,y  .weight; 
float  x_sum_total,y_sum_total; 
float  zero.limit; 

int  shade_low,shade_high; 


short  buffer[3100]; 
float  data[200000]; 
int  array.index  =  256; 
float  fft_data[256][256]; 
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off_t  seek_error,image_offset_bytes; 
unsigned  short  image_data[250][3200]; 
int  max_array_width  =  250; 
int  max_array_height  =  3000; 

int  win_buf; 

float  hist[260],hist_scale,max_hist; 
int  skip; 
int  buf_start; 

float  vscale; 
int  cxl,cx2; 
intcyl,cy2; 

int  fd,n_read,n_write, bytes;  /*  file  descriptors  */ 

FILE  *sfile,,"new_file,*old_file,*fp,*file_in,*fopen(),*fft_data_file; 

/+  ***************************************************************  */ 
int  isign; 


intnn[4]; 

j*  ***************************************************************  +/ 


/******«•  graphics  window  id’s  ***/ 
long  image_win,fftl_win,fft2_win,ana_win; 

main(argc.argv) 
int  argc; 
char  *argvD; 

{ 

#ifdef  show 
foregroundO; 

#endif 

r  =====================  REQUIRED  INPUTS  ======================== 

/*  A  35  mm  negitive  is  36mm  x  26  mm  ,  with  the  Nikon  scanner  at  a  pitch  of  1  */ 

/*  the  negative  is  digitized  to  125  pixels/mm  or  4500  x  3000  pixels  */ 

/*  For  use  on  the  Cray  the  size  of  the  buffer  and  picture  must  be  an  integer  */ 
/*  multiple  of  4  */ 

image_height  =  3000; 
image_width  =  4500; 
shadejow  =  10; 
shade_high  =  240; 

#ifdef  show 

printfl"  enter  Low  limit  of  the  picture  :  “); 
scanfr%d”,&shadeJow); 


printfT  Enter  the  High  limit  of  the  picture  :  “); 


scanfT%d”,&shade_high); 

#endif 

/*  The  program  must  know  the  diameter  of  a  particle  so  that  the  */ 

/*  zero  position  is  skipped  */ 
dia_part  =  6; 

/*  ===========================================================*/ 

/*  compute  the  number  of  windows  to  break  the  image  into.  This  will  be  the  number  */ 
/*  of  displacecent  vectors  per  image  ,  set  the  number  so  that  there  are  */ 

/*  about  8  particle  per  window  */ 

I*  At  this  time  the  value  must  be  either  16,32,64,128,or  256  *1 
xw_width  =  yw_width  =  64; 
nn[l]  =  128; 
nn[2]  =  128; 

ifl  nn[l]  >  array_index  t  I  nn[2]  >  array _index) 

( 

printfT  not  enough  memory  in  fit  arrays  \n”); 
exit(l); 

} 


iftimage  height  >  max_array_height ) 

{ 

printfT  height  too  big  \n"); 
exit(l); 

) 


iffxw.width  >  max_array  width  ) 

{ 

printfT  xw_width  too  big  \n"); 
exit(l); 

} 


/*  for  the  each  input  buffer  */ 
n_wind_sec_y  =  image_height/yw_width; 

n_wind_sec_x  =  image_width  /  xw_width;  /*  number  of  windows  in  the  x  direction  */ 
n_windows  =  n_wind_sec_x  *  n_wind_sec_y^*  total  number  of  windows  and  vectors  */ 

/***  open  the  Image  data  file  ***/ 
image_file  =  open(argv[l],0); 

/***  open  the  displacment  data  file  **/ 
new_file  =  fopen(argv[2],”w”); 

fprintfTnew_file,”%d  \n”,run_number); 
fprintflnew_file,"%d  \n”,frame_number); 
fprintflnew_file,"%d  \n”, image. width); 
fprintfTnew_file,”%d  \n”,image_h eight); 
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fprintfl[new_file,’’%d  \n",xw_ width); 
fprintf[new_file,”%d  \n”,yw_  width); 
fprintnnew_fi)e,”%d  \n”,n_wind_sec_x); 
fprintftnew_fx]e,”%d  \n”,n_wind_sec_y); 
fprintftnew_file,”%d  \n",nn[l]); 
fprintflnew_file,”%d  \n”,nn[2]); 

#ifdef  show 
vscale  =  5; 

cxl =  30; 

cx2  =  cxl  +  xw_width  *  vscale; 
cyl  =  30; 

cy2  =  cyl  +  yw_width  *  vscale; 

prefposition(cxl,cx2,cyl,cy2); 

image_win  =  winopen(“IMAGE  WINDOW  “); 

ortho2(0,xw_width,0,yw_width); 

RGBmodeO; 

gconfigO; 

cx2  =  1200  -  20; 

cxl  =  cx2  -  nn[l]  *  vscale; 

cy2  =  cyl  +  nn[2]  *  vscale; 

prefposition(cx  I,cx2,cyl,cy2); 
ffll_win  =  winopen(“  FFT  number  1”); 
ortho2(0,nn[l],0,nn[2]); 

RGBmodeO; 

gconfigO; 

cyl  =  cy2  +  30; 

cy2  =  cyl  +  nn[2]/2  *  vscale; 

prefposition(cxl,cx2,cyl,cy2); 
fft2_win  =  winopen(“  FFT  number  2”); 
ortho2(0,nn[l],0,nn[2]/2); 

RGBmodeO; 

gconfigO; 

#endif 


win_rel  =  1; 
x_offset  =  0; 

f*  compute  the  number  of  buffers  V 
buffer_bytes  =  xw_width  *  image_height  *  2; 

start_buffer  =  0; 

#ifdef  show 

printff"  enter  the  starting  buffer  :  “); 
scanfl“%d”,&start_buffer); 

#endif 
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end_buffer  =  n_wind_sec_x; 


/* 

start_buffer  =  30; 
endjauffer  =31; 

V 

for  (bufjoop  =  start_buffer;  bufjoop  <  endjiuffer;  bufJoop++) 

{ 

getJmageJileO; 

pivO; 

) 


/*  close  the  files  */ 

fclose(newjile); 

close(image_file); 


system(“  ps  -ef  I  grep  lapointe  »  acct.file  “); 

} 


get_image_file() 

[ 

printfl"  reading  data  An”); 

/*  Read  in  all  of  the  data  line  by  line  */ 

I*  the  data  is  read  into  the  array  “  buffer”  *1 

/*  Since  the  data  went  from  a  16  bit  machine  (PC)  to  a  32  bit  machine  (SGI)  */ 
/*  each  element  of  the  array  is  bit  shifted  by  8  bits  */ 

/*  The  data  on  the  PC  is  l’s  complement,  the  data  on  the  SGI  is  2’s  */ 

/*  complement,  so  if  the  data  is  less  than  0  add  256  */ 

image_offset_bytes  =  bufjoop  *  buffer  J>ytes; 

seek_error  =  lseek(image_file,image_offset_bytes,SEEK_SET); 
iftseek  error  ==  -1) 

( 

printff"  start  getting  data  seek  error  =  :  %ld  \n  “,seek_error); 
exit(l); 

) 


for  ( j  =  0  ;  j  <  xw_width  ;  ++j) 

{ 

n_read  =  read(imagejile, buffer, image Jieight  *  2); 
for(  k  =  0  ;  k  <  image Jieight;  ++k) 

( 

shade  =  buffer[k]  »  8; 
if(  shade  <  0  )  shade  +=  256; 
if(  shade  <  shade  Jo  w  ) 
shade  =  0; 

iff  shade  >  shadejiigh  ) 

shade  =  shadejiigh; 


G-8  • 


image_data[j][k]  =  shade; 


piv() 

{ 

ywl  =  0;  /*  y  */ 

win  =  1;  f*  first  window  number  */ 

/*  in  This  routine  loop  through  each  window,  get  the  two  comers  of  the  window  */ 

xwl  =0; 

xw2  =  xw_width; 

x_offset  =  buf_loop  *  xw_width  ; 

xwlout  =  x_offset; 

xw2out  =  xw2  +  x_offset; 

for  (y_win_step  =  0  ;  y  win  step  <  n_wind_sec_y;  y_win_step++) 

( 

ywl  =  y_win_step  *  yw_width; 
yw2  =  ywl  +  yw_width; 


fft_piv(); 

fprintflnew_file,”%d  %d  %d  %d  %d  %d  %d  %f  %f  %f  %f\n”, 

win_rel,xwlout,ywl,xw2out,yw2,peak_value_x,peak_value_y, 
avg_value_x  ,a  vg_value_y  ,fft_ratio  ,cen_ratio) ; 


win++; 

win_rel++; 

) 

) 

fft_piv() 

[ 

/**  Do  the  FFT  on  the  IMAGE  to  get  the  YOUNGS  FRINGES  */ 


tifdef  show 

show_image(); 

#endif 

pass_one(); 

two_d_fft(); 

fft_magnitude(); 

#ifdef  show 
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show_mag_l(); 

#endif 


/*  Do  the  FFT  on  the  FRINGES  to  get  the  DISPLACEMENT  */ 

pass_two(); 

two_d_fft(); 

fft_magnitude(); 

get_disp(); 

} 


/*  ========================================================= 

/*  ===================================================  v 

/*  =======================================================  */ 


pass_one() 

{ 

/*  On  the  first  pass  */ 

/*  Swap  the  2D  Real  data  input  array  to  a  1  D  */ 

/*  Complex  data  input  array  */ 

/*  zero  out  the  array  */ 
aindex  =  1; 

for(  k  =  0  ;  k  <  nn[2]  ;  k++) 

( 

for  (j  =  0  ;  j  <  nn[l];  j++) 

( 

data[aindex]  =  0; 
datalaindex  +  1  ]  =  0  ; 
aindex  +=  2; 

) 

} 

aindex  =  1; 

for(  k  =  0  ;  k  <  yw_width  ;  k++) 

( 

for  (j  =  0  ;  j  <  xw  width;  j++) 

{ 

datalaindex]  =  (float)image_data[j][k  +ywl]; 
datafaindex  +  1  ]  =  0  ; 
aindex  +=  2; 

} 

aindex  =  k  *  nn[l]  *2  +  1; 

} 

} 


/* 


/*  ================================================== 

pass_two() 

{ 

/*  On  the  second  Pass  ,  Take  the  Fringe  pattern  calculated  by  the  */ 
/*  first  2D  FFT  and  do  a  2DFFT  on  the  Fringe  Pattern  */ 

/*  This  will  give  a  displacment  vector  */ 


*/ 

*/ 

*/ 


aindex=  1 ; 


for(  k  =  0  ;  k  <  nn[2] ;  k++) 

( 

for  (j  =  0  ;  j  <  nn[l];  j++) 

( 

datalaindex]  =  fft_datafj][k]; 
data[aindex  +  1  ]  =  0  ; 
aindex  +=  2; 

) 


} 

/*=================================================== 

I*  =================================================== 

/*  =================================================== 

two_d_ffl() 

( 

float  theta_ratio; 
int  ntot,idim,ndim,n; 
int  il,i2,i3; 
int  ipl,ip2,ip3; 
int  ibit,i2rev; 
inti3rev,i3revi; 
intnrem,  nprev; 
intifpl,ifp2; 

int  kl,k2,k3; 

float  WT,wi,wpr,wpi,wtemp, theta; 
float  powx,  powy; 
float  tempi, tempr; 

/*  Set  up  the  sobroutine  for  a  2D  FFT  ,  ndim  =  2*1 
ndim  =  2; 
isign  =  1; 

/*  The  2D  FFT  is  taken  from  the  NUMERAICAL  RECIPES  BOOK  */ 
/*  The  codewas  converted  from  Fortran  to  C  */ 

ntot  =  1; 

for(  idim  =  1;  idim  <=  ndim;  idim++) 

( 

ntot  =  ntot  *  nn[idim]; 

} 

nprev  =  1; 

foK  idim  =  1;  idim  <=  ndim;  idim++) 

{ 

n  =  nn[idim]; 

nrem  =  ntot/(n  *  nprev); 

ipl  =  2  *  nprev; 

ip2  =  ipl  *  n; 

ip3  =  ip2  *  nrem; 

i2rev  =  1 ; 


*/ 

*/ 

*/ 
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for(  i2  =  1;  i2  <=  ip2;  i2  +=  ipl) 

{ 

if(  i2  <  i2rev) 

{ 

for(il  =  i2;  il  <=  i2+ipl  -2;  il  +=2) 

( 

for(i3  =  il;  i3  <=ip3;  i3  +=ip2) 

l 

i3rev  =  i2rev  +  i3  -  i2; 
tempr  =  data[i3]; 
tempi  =  data[i3  +  1]; 
data[i3]  =  data[i3rev]; 
data[i3  +1]  =  data[i3rev  +1]; 
data[i3rev)  =  tempr; 
data[i3rev  +1]  =  tempi; 

) 

) 

) 

ibit  =  ip2/2.0; 

while(  (ibit  >=  ipl)  &&  (i2rev  >  ibit)) 

( 

i2rev  =  i2rev  -  ibit; 
ibit  =  ibit/2.0; 

) 

i2rev  =  i2rev  +  ibit; 

) 

ifpl  =  ipl; 

while(  ifpl  <  ip2) 

{ 

ifp2  =  2.0  *  ifpl; 

theta  =  isign  *  6.283185  /  (ifp2  /  ipl); 

wpr  =  -2.0  *  sin(0.50  *  theta)  *  sin(0.50  *  theta); 

wpi  =  sin(theta); 

wr  =  1.0; 
wi  =  0.0; 

for(i3  =  1;  i3  <=  ifpl;  i3  +=ipl) 

( 

for(il  =  i3;  il  <=  i3  +  ipl  -2;  il  +=2) 

{ 

for  (i2  =  il;  i2  <=  ip3;  i2  +=  ifp2) 

{ 

kl  =  i2; 

k2  =  kl  +  ifpl; 

tempr  =  wt  *  data[k2]  -  wi  *  data[k2  +1]; 
tempi  =  wr  *  data[k2  +1]  +  wi  *  data[k2]; 
data[k2]  =  datafcl]  -  tempr; 
data[k2+l]  =  datatkl  +1]  -  tempi; 
data[kl]  =  data[kl]  +  tempr; 
data[kl+l]  =  datafkl+1]  +  tempi; 
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} 


) 

wtemp  =  wr; 

wr  =  wr  *  wpr  -  wi  *  wpi  +  wr; 
wi  =  wi  *  wpr  +  wtemp  *  wpi  +  wi; 
} 


ifpl  =  ifp2; 

) 

nprev  =  n  *  nprev; 

} 


) 


fft_magnitude() 

{ 


/*  Compute  the  Magnitude  of  the  FFT  */ 
aindex=  1  ; 

foK  k  =  0  ;  k  <  nn[2] ;  k++) 

{ 

for  (j  =  0  ;  j  <  nn[l];j++) 

{ 

square  =  datafaindex]  *  data[aindex]  +  data[aindex  +1  ]  *  data[aindex  +1]; 
magnitude  =  pow(square,0.50); 
fft_data|j][k]  =  magnitude  ; 
aindex  +=  2; 

} 


} 

/*==================================================================*/ 

/*  =========================================================  •/ 

I*  ========================================================  */ 


get_disp() 

( 

find_peak(); 

first_peak_value  =  fFt_peak; 
first_peak_x  =  peak_value_x; 
first_peak_y  =  peak_value_y; 


#ifdef  show 

show_mag_2(); 

#endif 


/*  zero  out  the  first  peak  */ 
j  =  peak_value_x; 
k  =  peak_value_y; 
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ajl  =  j  -  dia_part/2; 
iflajl  <  0)  ajl  =  0; 
aj2  =  j  +  dia_part/2; 
if[aj2  >  nn[l])  aj2  =  nn[l]; 

akl  =  k  -  dia_part/2; 
ifTakl  <  0)  akl  =  0; 
ak2  =  k  +  dia_part/2; 
iflak2  >  nn[l])  ak2  =  nn[l]; 

for(j  =  £yl;  j  <=  aj2;  j++) 

{ 

for(k  =  akl;  k  <=  ak2;  k++) 

( 

fft_data[j][k]  =  0; 

) 

) 


/*  find  the  second  peak  */ 
find_peak(); 

second_peak_value  =  fft_peak; 
second_peak_x  =  peak_value_x; 
second_peak_y  =  peak_value_y; 

fft_peak  =  first_peak_value; 
peak_value_x  =  first_peak_x; 
peak_value _y  =  first_peak_y; 


fft_ratio  =  first_peak_value  /  second_peak_value; 

iftpeak_value_x  >=  nn[l]/2) 

peak_value_x  =  peak_value_x  -  nn[l]; 

avg_value_x  =  peak_value_x; 
avg_value_y  =  peak_value_y; 


) 


V 
*/ 
*/ 

find_peak() 

{ 

peak_value_x  =  0; 
peak_value_y  =  0; 
fft_peak  =  0.0; 

/*  On  the  Second  Pass  Of  the  FFT  get  the  Peak  value  for  the  */ 

/*  Displacement  Vector  .  Since  the  2D  FFT  folds  over  ,  only  look  */ 

/*  at  half  of  the  map  */ 

for(  k  =  dia_part  ;  k  <  nn[2]/2  ;  k++) 


/* 

/* 
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{ 

for  (j  =  0  ;  j  <=  dia_part  ;  j++) 

{ 

iflfll_data[j][k]  >  fft_peak) 

{ 

fft_peak  =  fft_data[j]rk]; 
peak_value_x  = j; 
peak  value_y  =  k; 

] 


} 

) 

for(  k  =  0  ;  k  <  nn[2]/2  ;  k++) 

( 

for  (j  =  dia_part ;  j  <  nn[l]  -  dia_part  ;  j++) 

{ 

iflfft_data[j][k]  >  fft_peak) 

( 

fft_peak  =  fft_data[j][k]; 
peak_value_x  =  j; 
peak_value_y  =  k; 

} 

) 

} 

for(  k  =  dia_part  ;  k  <  nn[2]/2;  k++) 

( 

for  (j  =  nn[l]  -  dia_part ;  j  <=  nn[l]  ;  j++) 

( 

iflfift_data[j][k]  >  ffl_peak) 

( 

fft_peak  =  fft_data[j][k]; 
peak_value_x  =  j; 
peak  value_y  =  k; 

) 


} 


) 


) 


/* 

/* 

t* 

/* 

r 

i* 


*/ 

v 

v 

*/ 

V 

V 


#ifdef  show 

I*  ROUTINES  BELOW  THIS  LINE  ARE  ONLY  USED  TO  DEBUG  THE  CODE  **/ 
make_imageO 
( 
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int  al,a2,bl,b2; 
int  plx,ply,p2x,p2y,p3x,p3y; 
int  ptl; 
int  pt2; 

float  ratio_x,ratio_y; 
int  vel_x,vel  _y; 
int  ix,iy; 
int  npoints  ; 

run_n  umber  =  1234; 
frame_number  =  0; 
al  =  0; 

a2  =  xw_width/2; 
bl  =  0; 

b2  =  yw_width/2  ; 

printfT  a  %d  %d  b  %d  %d  \n",al,a2,bl,b2); 


/*  input  the  maximum  expected  displacment  in  the  x  and  y  directions  */ 
max_x_disp  =  40; 
max_y_disp  =  40; 

ratio_x  =  ((float)max_x_disp  -  103/1; 
ratio_y  =  ((f!oat)max_y_disp  - 10)/1; 

xwin  =  image_width; 
ywin  =  image_height; 

printfT  enter  x  vel :  “  ); 
scanfl“%d”,&vel_x); 
printfT  enter  y  vel :  “); 
scanft“%d”,&vel_y); 

printfT  number  of  points:  “); 
scanfT%d”,&npoints); 


fori  k  =  0  ;  k  <  npoints;  ++k) 

l 

ptl  =  (float)rand()/100; 
pt2  =  (float)rand(yiOO; 

while  (ptl  <  al  I  i  ptl  >  a2  I  I  pt2  <  bl  II  pt2  >  b2) 

{ 

ptl  =  (float)rand(yiOO; 
pt2  =  (float)rand(yiOO; 

] 

plx  =  ptl  +  j  *  xw_width; 
ply  =  pt2  +  i  *  yw_width; 
p2x  =  plx  +  vel_x  ; 
p2y  =  ply  +  vel_y  ; 
p3x  =  plx  +  vel_x  *2  ; 
p3y  =  ply  +  vel_y  *2  ; 
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Uplx; 
m  -  ply; 
make_dot(); 
1  =  p2x; 
m  =  p2y; 
make_dot(); 
1  =  p3x; 
m  =  p3y; 
make_dot(); 


) 


) 

make_dot() 

( 

image_data[l][m]  =  255; 

/* 

image_datatl-l][m]  =  200; 
image_data[l+l][m]  =  200; 
image_data[l][in+l]  =  200; 
image_data[l]ttn-l]  =  200; 

*/ 

} 

show  imageO 

( 

skip  =  1; 

winset(image_win); 
rgbvector[0]  =  200; 
rgbvector[l]  =  200; 
rgbvector[2]  =  0; 

printfT1  in  show  image  x  is  %d  %d  y  is  %d  %d  \n”,xwl,xw2,ywl,yw2); 

for  ( j  =  0  ;  j  <  xw_width  ;  j+=  skip  ) 

( 

foii  k  =  0  ;  k  <  yw_width;  k+=  skip) 

{ 

rgbvector[0]  =  image_data[j  +  xwl][k  +  ywl]; 
rgbvector[l]  =  image_data[j  +  xwl][k  +  ywl]; 
rgbvector[2]  =  image_data(j  +  xwl][k  +  ywl]; 

c3s(rgbvector); 

rectfi(j,kj+l,k+l); 

) 

} 


) 
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show_mag_l() 

( 

float  max_image; 
float  minjmage; 
float  slope,  b; 

minjmage  =  1000; 
max_image  =  0; 


winset(fftl_win); 

rgbvectortO]  =  200; 
rgbvectortl]  =  rgbvector[2]  =  0; 
c3s(rgbvector); 
clearQ; 


for  (j  =  0  ;j  <  nn[l]  ;j++  ) 

{ 

for(  k  =  0  ;  k  <  nn[2];  k++) 

{ 

ifT  fft_data[j][k]  >  maxjmage) 

{ 

max_image  =  fft_datafj][k]; 

) 

if!  fft_data[j][k]  <  minjmage) 

( 

minjmage  =  flft_data[j][k]; 

) 

) 

) 


slope  =  255.0  /  (maxjmage  -  minjmage); 
b  =  slope  *  minjmage; 

for  (j  =  0  ;j  <  nn[l];j++  ) 

{ 

for<  k  =  0  ;  k  <  nn[2];  k++) 

( 

cp  =  fft_data[j][k]  *  slope  -  b  ; 


/* 

get_color(); 

V 

rgbvectortO]  =  cp  ; 
rgbvectortl]  =  cp  ; 
rgbvector[2]  =  cp  ; 
c3s(rgbvector); 
rectfi(j,kj+l,k+l); 


G-18  ■ 


!* 

pnt2i(j,k); 

V 

) 


) 


show_mag_2() 

( 

float  max_image; 
float  min_image; 
float  slope,  b; 

min_image  =  1000; 
maxjmage  =  0; 

winset(fft2_win); 


rgbvector[0]  =  200; 
rgbvector[l]  =  rgbvector[2]  =  0; 
c3s(rgbvector); 
clearO; 

for  ( j  =  0  ;j  <  dia_part ;  j++  ) 

( 

for(  k  =  dia_part ;  k  <  nn[2]/2;  k++) 

{ 

ift  fft_data[j][k]  >  max  image) 

( 

maxjmage  -  fft_data|j][k]; 

1 

iR  fft  data[j][k]  <  min  image) 

( 

min  Jmage  =  fft_data[j][k]; 


) 

for  ( j  =  dia_part ;  j  <  nn[l]  -  dia_part ;  j++  ) 

( 

for(  k  =  0  ;  k  <  nn[2]/2;  k++) 

( 

ifl  fH_data|j][k]  >  maxjmage) 

( 

maxjmage  =  fft_datalj][k]; 

} 

iR  fft_data[j][k]  <  min  Jmage) 

( 
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minjmage  =  flft_data[j][k]; 
] 


/* 

printfT  max  image  =  %f  min  =  %f  \n”, maxjmage,  minjmage); 

maxjmage  =  2000000.0; 
min  image  =  10000.0; 

*/ 

slope  =  255.0  /  (maxjmage  -  minjmage); 
b  =  slope  *  minjmage; 


for  ( j  =  0  ;  j  <  nn[l] ;  j++  ) 

( 

for(  k  =  0  ;  k  <  nn[2];  k++) 

{ 

cp  =  fft_data[j][k]  *  slope  -  b; 

get_color(); 

c3s(rgbvector); 


rectfi(j  ,k  j + l,k+ 1) ; 
f* 

*/ 


pnt2i(j,k); 

1 


) 

get_color<) 

{ 


rgbvector[l]  =  cp  ; 
rgbvector[l]  =  0; 
rgbvector[2]  =  0; 

iftcp  <  250) 

{ 

rgbvector[0]  =  0; 
rgbvector[l]  =  cp  *  255.0/250.0  ; 
rgbvector[2]  =  0; 

1 

iflcp  <  245) 

{ 

rgbvector[0]  =  0; 
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rgbvector[l]  =  0; 
rgbvector[21  =  cp  *  255.0/245.0; 

} 

ifTcp  <  240) 

{ 

rgbvectorfO]  =  0; 
rgbvector[l]  =  cp  *  255.0/240.0; 
rgbvectorfO]  =  cp  *  255.0/240.0  ; 
) 

ifTcp  <  230) 

{ 

rgbvectorfO]  =  cp  *  255.0/230.0; 
rgbvectorfl]  =  0; 
rgbvector[2]  =  cp  *  255.0/230.0  ; 
) 

ifTcp  <  220) 

{ 

rgbvectorfO]  =  cp ; 
rgbvectorfl]  =  cp ; 
rgbvector[2]  =  cp  ; 

) 

ifTcp  <  10) 

l 

rgbvectorfO]  =  0 ; 
rgbvectorfl]  =  0  ; 
rgbvector[2]  =  0  ; 

] 


) 

#endif 
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Appendix  H 


2D  FFT  PIV  CODE  FOR  THE  CRAY:  CRAYJFFTPIV.C 


H-l/H-2 
Reverse  Blank 


/**  PARTICLE  IMAGE  VELOCEMETRY  PROGRAM  **/ 
/+*  Using  the  2D  fit  method  **/ 

/**  By :  Kenneth  LaPointe  **/ 

/**  NUSC,  Code  8322  **/ 

/**  Last  rev  as  of  Aug  17  1992  **/ 

/**  This  program  reads  in  data  files  and  image  **/ 

/**  files  to  use  an  2D  FFT  program  **/ 

/**  to  compute  the  displacemnet  between  **/ 

/**  particles  **/ 

/**  compile  on  CRAY  ONLY  */ 

/**  Run  this  Script  to  compile  the  code  on  the  Cray 
and  move  it  to  the  working  directory 

echo  Copy  the  Source  code  to  the  CRAY 
rep  fltpiv.c  cray:fitpiv.c 

echo  compile  The  source  code 
rsh  cray  cc  -03  -c  fftpiv.c 

echo  Link  the  code 

rsh  cray  segldr  -lm  fftpiv.o  -o  fftpiv 

echo  Make  the  tmp  user  directory 
rsh  cray  mkdir  /usr/tmp/PFV 

echo  move  the  executable  code  to  the  tmp  dir 
rsh  cray  mv  fftpiv  /usr/tmp/PIV 

********************  */ 

#include  <stdio.h> 

#in elude  <math.h> 

#include  <unistd.h> 

#include  <sy&/types.h> 

#include  <complex.h> 

int  end_buffer; 

int  fff_loop  ; 

int  fft_pass; 

intij,k,l,m,n,o,p; 

int  aj  I,aj2,aj3,ak  I,ak2,ak3; 

float  a,b,c,d.e,f,g; 

double  square,  magnitude; 

float  fft_peak; 
float  fff_ratio; 
float  cen_ratio; 


float  firstjpeak_value; 
int  first_peak_x; 
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int  first_peak_y; 

float  second_peak_value; 
int  second_peak_x; 
int  second_peak_y; 
int  aindex; 

int  junk; 

int  nbuff; 

int  buffer_bytes; 

int  buf_loop,x_offset,start_buffer; 
int  pic_width; 
int  win_rel; 

short  rgbvector[3]; 
long  datastore; 

int  xwl,ywl;  /*  lower  left  corner  of  2  D  array  is  xwl.ywl  */ 

int  xw2,yw2;  /*  upper  right  comer  of  2  D  array  is  xw2,yw2  */ 

int  xw_width,yw_width; 

int  n_windows,win; 
int  max_frames,run_number; 
int  x_win_step,y_win_step; 
int  total; 

int  xwlout,xw2out; 

int  max_x_disp,max_y_disp; 

int  image_file; 

int  image_height,image_width; 

char  old_disp_name[50]; 
char  new_disp_name[50]; 

int  shade; 

int  n_wind_sec_x,  n_wind_sec_y; 

int  frame_number; 

int  dia_part; 

intxwin,ywin; 

float  skipx,  skipy; 


int  draw_data; 
int  shift; 

int  peak_value_x; 
int  peak_value_y; 
float  avg_value_x; 
float  avg_value_y; 

float  x_sum  ,y_sum  ,x_weight,y_weight; 
float  x_sum_total,y_sum_total; 
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float  zero_limit; 

int  shade_low,shade_high; 

longbuffer[1024]; 

unsigned  long  sl,s2,s3,s4; 

int  array_index  =  300; 

complex  double  complex_data[300][300]; 

float  fft_data[300][300]; 

float  table[2048]; 

int  n  table; 

off_t  seek_error,image_offset_bytes; 
unsigned  short  image_data[256][3200]; 
int  max_array_width  =  256; 
int  max_array_h  eight  =  3000; 

int  win_buf; 

int  skip; 
int  buf_start; 

int  fd,n_read,n_write, bytes;  /*  file  descriptors  */ 

FILE  *new_file,*old_file,*fopen(); 

int  isign; 

int  n_fft_x,n_fft_y; 


main(argc,argv) 
int  argc; 
char  *argv[]; 

( 

printft"  Start  of  Program  \n”); 

/*  =====================  REQUIRED  INPUTS  ============================  */ 

/*  A  35  mm  negitive  is  36mm  x  26  mm  ,  with  the  Nikon  scanner  at  a  pitch  of  1  */ 

/*  the  negative  is  digitized  to  125  pixels/mm  or  4500  x  3000  pixels  */ 

/*  For  use  on  the  Cray  the  size  of  the  buffer  and  picture  must  be  an  integer  */ 

/*  multiple  of  4  */ 

image_height  =  3000; 
image_width  =  4500; 
shadejow  =  10; 
shade_high  =  240; 

f*  The  program  must  know  the  diameter  of  a  particle  so  that  the  */ 

/*  zero  position  is  skipped  */ 
dia_part  =  6; 

/*  ===========================================================*/ 

/*  compute  the  number  of  windows  to  break  the  image  into.  This  will  be  the  number  */ 

/*  of  displacecent  vectors  per  image  ,  set  the  number  so  that  there  are  */ 

/*  about  8  particle  per  window  */ 
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xw_width  =  yw_width  =  40; 


/*  The  number  of  points  in  the  ffl  value  should  be  either 
I*  16, 32,64,128, or  256 

/*  the  value  must  be  greater  than  the  window  width 
n_fft_x  =  128; 
n_fft_y  =  128; 

/* 

printfT  Check  memory  allocation  \n”); 

ift  n_ffl_x  >  array  index  I  I  n_fft_y  >  array _index) 

{ 

printf(u  not  enough  memory  in  fft  arrays  \n”); 
exit(l); 

} 


iffimage_height  >  max_array_height ) 

( 

printfl"  height  too  big  \n”); 
exit(l); 

) 

iffxw_width  >  max_array  width  ) 

{ 

printfll"  xw_width  too  big  \n”); 
exit(l); 

) 

*/ 


f*  for  the  each  input  buffer  */ 
n_wind_sec_y  =  image_height/yw_width; 

n_wind_sec_x  =  image_width  /  xw_ width;  /*  number  of  windows  in  the  x  direction  */ 
n_windows  =  n_wind_sec_x  *  n_wind_sec_y;  /*  total  number  of  windows  and  vectors  */ 

printft"  Open  Data  Files  \n”); 

/***  open  the  Image  data  file  ***/ 
image_file  =  open(argv[l],0); 

/***  open  the  displacment  data  file  **/ 
new_file  =  fopen(argv[2],”w”); 

fprintffne w_file,”%d  \  n”,run_n umber); 
fprintfinew_file,”%d  \n”,frame_number); 
fprintffnew_file,”%d  \n”,image_width); 
fprintfinew_file,”%d  \n”,image_h eight); 
fprintffnew_file,”%d  \n”,xw_width); 
fprintftnew_file,”%d  \n”,yw_width); 
fprintftnew_file,”%d  \n",n_wind_sec_x); 
fprintfl[new_file,”%d  \n”,n_wind_sec _y); 


*/ 

V 

*/ 
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fprintf(new_file,”%d  \n”,n_flft_x); 
fprintftnew_file,”%d  \n”,n_fft_y); 

win_rel  =  1; 
x_oflfset  =  0; 

/*  compute  the  number  of  buffers  */ 
buffer_bytes  =  xw_width  *  image_height  *  2; 


printfT  Start  The  x  loop  \n”); 
start_buflfer  =  0; 
end_buffer  =  n_wind_sec_x; 

for  (buf_loop  =  start_buffer;  buf  loop  <  endjouflfer;  buf_loop++) 

( 

get_image_file(); 

piv(); 

1 


/*  close  the  files  */ 

fclose(new_file); 

close(image_file); 


system(“  ps  -ef  I  grep  lapointe  »  acct.file  “); 

1 


get_image_file() 

( 

/*  Read  in  all  of  the  data  line  by  line  */ 

/*  the  data  is  read  into  the  array  “  buffer”  */ 

/*  Since  the  data  went  from  a  16  bit  machine  (PC)  to  a  64  bit  machine  (CRAY)  */ 
/*  1  word  on  the  cray  contains  4  PC  words  */ 

image_oflfset_bytes  =  bufjoop  *  buffer_bytes; 

seek_error  =  lseek(image_fUe,image_oflfset_bytes,SEEK_SET); 

for  ( j  =  0  ;  j  <  xw_width  ;  ++j) 

{ 

m  =  0; 

n_read  =  read(image_file,bufFer,image_height  *  2); 
foK  k  =  0  ;  k  <  image_height/4;  k++) 

{ 

image_data[j][m]  =  image_data(j][m+l]  =  image_data(j][m+2]  = 
image_data[j][m+3]  =  0; 
si  =  bufferfk]  ; 
image_data[j][m  ]  =  si  »  56; 

iflimage_data[j][m]  <  shade_low  I  I  image_data(j][m]  >  255) 
image_data[j][m]  =  0; 

s2  =  buffer[k]  «  16; 
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image_data[j][in+l]  =  s2  »  56; 

iflimage_data[j][in+l]  <  shadejow  I  I  image_data[j][m+l]  >  255) 
image_data(j][m+l]  =  0; 

s3  =  bufTer[k]  «  32; 

image_data[j][m+2]  =  s3  »  56; 

ifIimage_data(j][m+2]  <  shadejow  I  I  image_data[j][m+2]  >  255) 
image_data[j][m+2)  =  0; 

s4  =  bufferfk]  «  48; 

image_data[j][m+3]  =  s4  »  56; 

ifTimage_data[j][m+3]  <  shadejow  I  I  image_dataij][m+3]  >  255) 
image_data[j][m+3]  =  0; 


m  +=4; 

) 

) 

) 


piv() 

( 

ywl  =  0; 

win  =  1; 

xwl  =  0; 

xw2  =  xw_ width; 

x_offset  =  bufjoop  *  xw_width  ; 

xwlout  =  x_offset; 

xw2out  =  xw2  +  x_offset; 

for  (y_win_step  =  0  ;  y_win_step  <  n_wind  sec_y;  y_win_step++) 

( 

ywl  =  y_win_step  *  jrw_width; 
yw2  =  ywl  +  yw_width; 

ffl_piv(); 

fprintfTnew_file,”%d  %d  %d  %d  %d  %d  %d  %f  %f  %f  %f  \n”, 

win_rel, xwlout, ywl, xw2out,yw2,peak_value_x,peak_value_y, 
avg_value_x,avg_value_y,ffl_ratio,cen_ratio); 

win++; 

win_rel++; 

) 

) 

fft_piv() 

( 

/**  Do  the  FFT  on  the  IMAGE  to  get  the  YOUNGS  FRINGES  */ 

pass_one(); 

two_d_fft(); 

fft_magnitude(); 
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I*  Do  the  FFT  on  the  FRINGES  to  get  the  DISPLACEMENT  */ 

pass_two(); 

two_d_fft(); 

fft_magnitude(); 

get_disp(); 


/*  ===================================== 

/*  ===================================== 

/*  ===================================== 

pass  oneO 

{ 

/*  On  the  first  pass  set  the  array  to  zero  */ 
for(  k  =  0  ;  k  <  n_fft_y  ;  k++) 

( 

for  (j  =  0  ;  j  <  n_fft_x;j++) 

{ 

complex_data[j][k]  =  0  +  Oi; 

) 

J 

/*  Then  input  the  data  into  the  complex  array  */ 
for ( k  =  0  ;  k  <  yw  width  ;  k++) 

( 

for  (j  =  0  ;j  <  xw_width; j++) 

{ 

complex_datafk][j]  =  (float)image_data[j][k  +  ywl]  +  Oi; 

) 

) 


*/ 

*/ 

*/ 


) 

/*================================================== 

/*  ================================================== 

/*================================================== 

pass_two() 

( 

/*  On  the  second  Pass  ,  Take  the  Fringe  pattern  calculated  by  the  *t 
/*  first  2D  FFT  and  put  it  into  the  input  array  for  the  next  pass*/ 


fort  k  =  0  ;  k  <  n  fft_y  ;  k++) 

{ 

for  (j  =  0  ;  j  <  n_flft_x;  j++) 

{ 

complex_data[k][j]  =  fft_data[j][k]  +  Oi; 

) 

) 

) 


*/ 

*/ 

*/ 
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/*  ==========================================================  */ 

two_d_fft() 

( 

float  work[263000]; 
float  scale; 

int  isign,nl,n2; 
int incrl; 
int  nwork; 
int  incr; 


isign  =  1; 
nl  =  n_fft_x; 
n2  =  n_flft_y; 
scale  =  1.0; 
incr  -  1; 
skip  =  n_f¥t_x; 

ntable  =  2048; 


nwork  =  262144; 


/*  See  the  CRAY  C  manaual  for  the  discussion  of  this  subroutine  */ 
CFFT2D(&isign,&nl,&n2,&scale,complex_data,&incr,&array_index, 

complex_data,&incr,&array_index, table, &ntable, work, &nwork); 


/* 

/* 


/* 


fft_magnitude() 

( 


/*  Compute  the  Magnitude  of  the  FFT  */ 
for(  k  =  0  ;  k  <  n_fft_y  ;  k++) 

( 

for  (j  =  0  ;  j  <  n_fft_x;  j++) 

{ 

magnitude  =  cabs(complex_data[k][j]); 
fft_data[j][k]  =  magnitude ; 

) 

) 


*/ 

V 

*/ 


) 


/*  =========================================================  V 


get_disp() 

{ 

/*  find  highest  peak  in  the  2D  FFT  data  V 
find_peak(); 
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first_peak_value  =  fft_peak; 
first_peak_x  =  peak_value_x; 
first_peak_y  =  peak_value_y; 


/*  zero  out  the  first  peak  data  to  find  the  next  highest  peak  */ 
j  =  peak_value_x; 
k  =  peak_value_y; 


ajl  =  j  -  dia_part/2; 

iftajl  <  0)  ajl  =  0; 

aj2  =  j  +  dia_part/2; 

ifiaj2  >  n_fFt_x)  aj2  =  n_fift_x; 

akl  =  k  -  dia_part/2; 
ifiakl  <  0)  akl  =  0; 
ak2  =  k  +  dia_part/2; 
i(Iak2  >  n_fft_y)  ak2  =  n_flt_y; 

for(j  =  ajl;  j  <=  aj2;  j++) 

( 

for(k  =  akl;  k  <=  ak2;  k++) 

{ 

fft  data[j][k]  =  0; 

) 

} 


/*  find  the  second  highest  peak  in  the  2D  FFT  data*/ 
find_peak(); 

second_peak_value  =  fft_peak; 
second_peak_x  =  peak_value_x; 
second_peak_y  =  peak_value_y; 


fft_peak  =  first_peak_value; 
peak_value_x  =  first_peak_x; 
peak_value_y  =  first_peak_y; 


fFt_ratio  =  first_peak_value  /  second_peak_value; 

/*  Decide  weather  the  peak  is  in  the  first  or  second  quadrant  */ 
ifTpeak_value_x  >=  n_fft_x  /  2) 

peak_value_x  =  peak_value_x  -  n_flft_x; 

/*  For  now  the  avg_peak  is  the  same  as  the  peak  value  */ 

/*  Later  code  will  be  added  to  compute  the  average  displacement  around  the  peak  */ 
avg_value_x  =  peak_value_x; 
avg_value_y  =  peak_value_y; 


} 
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/* ================================================== 

find_peak() 

( 

peak_value_x  =  0; 
peak_value_y  =  0; 
fft_peak  =  0.001; 

/*  After  the  Second  Pass  of  the  2D  FFT  get  the  Peak  value  for  the  */ 
/*  Displacement  Vector  .  Since  the  2D  FFT  folds  over  ,  only  look  */ 

/*  at  half  of  the  map  */ 

for(  k  =  0  ;  k  <  dia_part ;  k++) 

( 

for  (j  =  dia_part ;  j  <  n_fft  x  -  dia_part  -1  ;  j++) 

( 

iflfft_data[j][k]  >  fft_peak) 

{ 

fft_peak  =  fft_data[j][k]; 
peak_value_x  =  j; 
peak_value_y  =  k; 

) 


} 

} 


for(  k  =  dia_part ;  k  <  n  fft_y/2  -1;  k++) 

( 

for  (j  =  0;  j  <  n_fft_x  -1  ;  j++) 

{ 

ifTfft_data[j][k]  >  fft_peak) 

{ 

fft_peak  =  fft_datalj][k]; 
peak_value_x  =  j; 
peak_value_y  =  k; 

) 


Appendix  I 


VECTOR  DISPLAY  CODE:  VECTOR.C 


I-1/I-2 

Reverse  Blank 


/**  PARTICLE  IMAGE  VELOCEMETRY  DISPLAY  PROGRAM 
/**  Draw  vectors  on  screen  and  in  file  for  hp  plotter  **/ 

/**  By :  Kenneth  LaPointe  **/ 

/**  NUSC,  Code  8322  **/ 

/**  REV  C.  ,  Last  rev  as  of  20  Feb  1992  **/ 

/**  to  compile  type  :  cc  vector.c  -lgl_s  -lm  -o  vector  **/ 

#include  <gl/gl.h> 

#include  <stdio.h> 
tinclude  <math.h> 

int  true  =  1; 

int  false  =  0; 

int  window_too_small; 

int  ij,k,l,m,n,o,p; 

int  junk; 

int  win_width_x; 

int  win_width_y; 

int  xw_width,yw_width; 

int  n_windows,window_number; 

int  frame,  max_frames,run_number,frame_step; 

int  x_win_step,y_win_step,win_steps; 

int  total,max_corr; 

int  old_x_peak,  old_y_peak; 

int  xd_peak,yd_peak; 

int  disp_shift; 

float  avg_xd,avg_yd; 

float  x_end,y_end,x_center,y_center; 

float  magnitude; 

float  angle  ,  yx,  h_angle; 

float  xpl,xp2,xp3,xp4,ypl,yp2,yp3,yp4; 

int  winx.winy; 

float  a,b,c,d,e,f; 

float  rl,r2,r3,T4,r5,r6,r7,r8; 

float  rll.rhl; 

int  dal,da2,da3,da4; 

float  window_x,window_y; 
float  fft_ratio; 
float  grid_ratio; 

char  file_name[30];  /*  file  name  character  storage  *1 

char  title[30]; 

char  last_data_file[20]; 

int  image_height,image_width; 
int  xoffset; 


char  input_name[100]; 


int  shade; 

int  n_wind_sec_x,  n_wind_sec_y; 
int  frame_number; 
int  disp_file; 
int  lstep; 

int  xwl,xw2,ywl,yw2; 

int  fd,n_read,n_write, bytes;  /*  file  descriptors  */ 

float  win_scale; 

int  cmap; 

int  red,green,blue; 

int  wn; 

int  x_peak,y_peak; 

float  vector  [15000][6];  /*  0  magnitude 

1  angle 

2  x  center 

3  y  center 

4  fft_ratio 

5  quality  */ 

int  xfft,yfft; 
float  max_disp; 
int  ix; 

float  scale_factor; 
float  dim_ratio; 
float  laser_fTeq ; 
int  xmar,ymar; 

float  max_mag,  mratio.ncolors; 
int  pens,pratio; 
int  pen_number; 
int  xa,xb; 

int  showflag ,  closeflag.ratioflag; 
int  part_dia; 

float  mO,aO,inv,av,lowl,highl,error_limit; 
float  ratio Jimit; 

FILE  *sfile,*dsp_file,*old_fi]el*fp,*file_in,*hpI*fopen(); 

j*  ********t**************4'***+**#****  +  +  *****'t'**+**********>li+****+  */ 


main(argc.argv) 
int  argc; 
char  *argv[]; 

( 

foregroundO; 

filejn  =  fopen(“holder.dat”  ,”w"); 

fprintfTfile_in,”%s”,argv[l]); 

fclose(file_in); 

file_in  =  fopen(“holder.dat”  ,"r”); 


fscan  fXfile_in,  "%s”,file_name); 
fclose(file_in); 

disp_shift  =  0; 

cmap  =  1000; 

ncolors  =  10; 

pens  =  6; 

scale_factor  =  4; 

h.angle  =  5.  *  (2.  *  3.14  /  360.0); 

max_disp  =  100.0; 


printPT  enter  the  ratio  limit :  (  0  shows  all  vectors)  “); 
scanfl[“%r,&ratio_limit); 

part_dia  =  -1; 
ift  ratiojimit  !=  0) 

( 

printf(“  enter  the  particle  diameter  (pixels) :  “); 
scanf(“%d”,&part_dia); 

) 

/*  set  up  the  limit  for  neighboring  vectors  */ 
error_limit  =  0.25; 
lowl  =  1.0  -  error_limit; 
highl  =  1.0  +  errorjimit; 


/*  set  the  ratio  to  get  from  pixels  to  fVsec  */ 
dim_ratio  =  0.0520  ;  /*  ft/sec  per  pixel  */ 


printfT  OPEN  DATA  FILE  HEADER  %s  \n”, filename); 
dsp_file  =  fopen(file_name,V); 
read_datafile_header(); 
setup_plotter(); 


printf(“  Read  in  the  DATA  FILE  \n"); 

xwl  =  0; 
ywl  =  0; 
xoffset  =  0; 

for(wn  =  0  ;  wn  <  n_windows;  wn++) 

{ 

get_disp(); 

) 

fclose(dsp_file); 


/*  Set  up  the  color  scale  */ 

mratio  =  ncolors  /  (max_mag  *  scale_factor); 

pratio  =  pens  /  (max_mag  *  scale_factor); 
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printfT  colors  %f,  max  %f  ratio  %f  \n”,ncolors,max_mag,inratio); 


select_data(); 

interpolateO; 


open_vector_window(); 

iflratio  limit  ==  0) 

{ 

for(wn  =  0  ;  wn  <  n_windows;  wn++) 

( 

color(l); 

compu  te_vector() ; 

screen_arrow(); 

paper_arrow(); 


} 

else 

( 

foKwn  =  0;  wn  <  n_windows;  wn++) 

( 

if!vector[wn][5]  ==1) 

{ 

color(O); 

c  ompute_vector() ; 

screen_arrow(); 

paper_arrow(); 

1 

if(vector[wn][5]  ==  2) 

( 

color(l); 

compute_vector(); 

screen_arrow(); 

paper_arrow(); 

) 

} 

) 

fclose(hp); 


printfT"  finished  dumping  plot  data  \n”); 
sleep(960); 

) 


read_datafile_header() 
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( 


fscanftdsp_file,"%d  “,&run_number); 
fscanfldsp_fde,”%d  “,&frame_number); 
fscanfldsp_fde,”%d  “,&image_width); 
fscanfldsp_file,”%d  “,&image_height); 
fscanfldsp_file,”%d  “,&xw_width); 
fscanfldsp_file,”%d  “,&yw_width); 
fscanfldsp_file,”%d  “,&n_wind_sec_x); 
fscanfldsp_file,”%d  “,&n_wind_sec_y); 
fscanfldsp_file,”%d  “,&xfft); 
fscanfldsp_file ,”%d  “,&yfft); 

n_windows  =  n_wind_sec_x  *  n_wind_sec_y; 
ix  =  n_wind_sec_y; 


/* 

printfl"  %d  \n”,run_number); 

printf(“  frame  number  %d  \n",frame_number); 

printfl"  image  width  %d  \n”,image_width); 

printfl"  image  height  %d  \n",image_height); 

prints"  number  of  window  secs  in  x  %d  \n",n_wind_sec_x) 

printfl"  in  y  %d  \n”,n_wind_sec_y); 

*/ 

} 


setup_plotteT() 

{ 

/*  Open  the  HP  plotter  file  */ 
printfl"  open  the  Plotter  file  and  set  it  up  \n”); 
hp  =  fopen(“vector.hp”,”w”); 
fprintflhp,”SPl;  \n”); 
fprintflhp ,”PT.3;  \n”); 
fprintflhp, ”DT\03;  \n”); 


winx  =  image_width  *  11/8.5; 
winy  =  image_h  eight; 

lstep  =  0.02  *  winx; 

a  =  1.0/125.0; 

fprintflhp, ”SC  -100,%d,-100,%d;  \n”,winx,winy); 

/*  Draw  Run  ID  info  */ 
xmar  =  winx; 

ymar  =  image_height  - 100; 

fprintflhp,”  PU%d,%d;DI0,-l;LBFILE  NAME  :  %s  \03  \n”, 
xmar  ,ymar,file_name); 

fprintflhp,”  PU%d,%d;DI0,-l;LBRUN  ID  NUMBER  :  %d  \03  \n”, 
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xmar  -  lstep,ymar,run_number); 


fprintfThp,”  PU%d,%d;DIO,-l;LBFrame  Number  :  %d  \03  \n”, 
xmar  -  2  *  lstep,ymar,frame_number); 

fprintfThp,”  PU%d,%d;DIO,-l;LBPicture  size  :  %d  X  %d  pixels  \03  \n”, 
xmar  -  3  *  lstep,ymar,image_width,image_height); 

fprintfThp,”  PU%d,%d;DIO,-l;LBWindow  size  :  %d  X  %d  pixels  \03  \n”, 
xmar  -  4  *  lstep,ymar,xw_width,yw_width); 

1  =  image_width  /  xw_width; 
m  =  image_height/yw_width; 
n  =  1  *  m  ; 

fprintfThp,"  PU%d,%d;DI0,-l;LBVector  Count  :  %d  X  %d  :  total  %d  \03  \n”, 
xmar  -  5  *  lstep,ymar,l,m,n); 

fprintfThp,”  PU%d,%d;DIO,-l;LBCon  version  :  1  pixel  =  %f  ft/sec  \03  \n”, 
xmar  -  6  *  lstep,ymar,dim_ratio); 

fprintfThp,”  PU%d,%d;DIl,0  ;LB  +X  axis  \03  \n”,image_width  +  lstep.lstep); 
fprintfThp,”  PU%d,%d;DIO,-l;LB  +Y  axis  \03  \n”,-lstep,image_height  -  lstep); 

/*  draw  box  */ 

fprintfThp,”  PU%d,%d;  \n",0,0); 

fprintfThp,”  PD%d,%d;  \n",image_width  +  5  *  lstep, 0); 

fprintfThp,”  PU%d,%d;  \n”,0,0); 

fprintfThp,”  PD%d,%d;  \n”,0,image_height  +  5  *  lstep); 
fprintfThp,”  PU%d,%d;  \n”,image_width,0); 
fprintfThp,”  PD%d,%d;  Vn”,image_width,image_height); 
fprintfThp,”  PD%d,%d;  \n”,0,image_height); 

fprintfThp,”  PU%d,%d;DI,0,-l;LBScale:  50  pixels  =  %ffVsec  \03  \n”, 
image_width  +  lstep,image_height  -  4  *  lstep,  50.0  *  dim_ratio); 

fprintfThp,”PT.l;  \n”); 

xpl  =  image_width  +  lstep; 
ypl  =  image_height  -  2  *  lstep; 

wn  =  n_windows  +  1; 

I*  make  the  vector  scale  key  */ 
vector[wn][0]  =  50  *  scale_factor; 
vector[wn][l]  =  0; 

vector[wn][2]  =  image_width  +  lstep/2.; 
vector[wn][3]  =  image_height  -  2  *  lstep; 
compute_vector(); 
paper_arrow(); 

vector[wn][l]  =  90  *  2  *  3.14159  /  360; 
com  pute_vector() ; 
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paper_arrow(); 


} 

open_vector_window() 

( 

printfT  FLOW  VISUALIZATION  PROJECT  ,  DISPLAY  PROGRAM  \n”); 

prefposition(10,4500/3.5  ,10,3000/3.5); 

winopen(“  TEST  IMAGE  “); 

make_cmapO; 

ortho2(0,4500, 0,3000); 

color(7); 

clearO; 

color(3); 

} 

get_dispO 

( 


fscanfTdsp_file,”%d  %d  %d  %d  %d  %d  %d  %f  %f  %f  %f  \n  “,&window_number, 

&xwl,&ywl,&xw2,&;yw2,&x_peak,&y_peak,&avg_xd,&;avg_yd,&fFl_ratio,&grid_ratio); 

/*  use  the  peak  values  */ 
avg_xd  = x_peak; 
avg_yd  = y_peak; 

/*  get  rid  of  the  particle  size  limited  magnitudes  */ 

iflavg_xd  ==  part_dia  I  I  avg_yd  ==  part_dia) 

{ 

avg_xd  =  0; 
avg_yd  =  0; 

) 


magnitude  =  pow(avg_xd  *  avg_xd  +  avg_yd  *  avg_yd,.5); 

/*  Assume  a  +x  and  +y  flow  */ 
iflavg_xd  ==  0 ) 

( 

angle  =  1.57;  l*  90  deg  up  */ 

) 

else 

{ 

yx  =  avg_yd/avg_xd; 
angle  =  fatan(yx); 

) 

ifl  angle  <=  0) 

angle  =  angle  +  3.14159; 


vector!  wn]fO]  =  magnitude  *  scale_factor; 
vectorlwn][l]  =  angle; 
vector[wn][2]  =  (xwl  +  xw2)/2.0; 
vector[wn]f3]  =  (ywl  +  yw2)  /2.0; 
vector[wn][4]  =  fftjratio; 
vector[wn][5]  =  0; 


/*  get  the  peak  magnitude  for  the  color  scale  */ 
if  (magnitude  >  max_mag) 

max.mag  =  magnitude; 


select_data() 

{ 

for  (wn  =  0;  wn  <  n_windows;  wn++) 

( 

showflag  =  closeflag  =  ratioflag  =  0; 

/*  test  for  the  peak  to  peak  ratio  */ 
if!vector[wn][4]  >=  ratio Jimit ) 
ratioflag  =  100; 

/*  Test  for  a  vectors  with  the  close  to  same  magnitude  and  angle  V 
mO  =  vector[wn][0]; 
aO  =  vector[wn][l]; 

/**  Loop  through  the  3x3  matrix  of  vectors  centered  around  vector! wn]  */ 
for(j  =  -1 ;  j  <=  1 ;  j++) 

( 

for(k  =  -ix  ;  k  <=  ix  ;  k  +=  ix) 

( 

if(  j  +  k  !=  0) 

( 

mv  =  vector!  wn  +  j  +  k][0]; 
av  =  vector!  wn  +  j  +  kit  1]; 

ifl  ((mv  <=  mO  *  highl)  &&  (mv  >=  mO  *  lowl))  && 
((av  <=  aO  *  highl)  &&  (av  >=  aO  *  lowl))  ) 

( 

closeflag  +=  3; 

} 

) 

) 

) 

showflag  =  ratioflag  +  closeflag; 

/*  greater  than  confidence  ratio  and  more  than  1  close  neighbors  */ 
if( showflag  >  103) 

( 
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vector[wnl[5]  =  1; 

) 


} 


) 


interpolateO 

( 

float  sum_mag; 
float  sum_ang; 
float  sum_vec; 


/*  Interpolate  where  there  is  no  data  */ 
for(wn  =  0  ;  wn  <  n_windows;  wn++) 

( 


sum_mag  =  sum_ang  =  sum_vec  =  0; 
iflvector[wn][5]  ==  0) 

{ 

/**  Loop  through  the  3x3  matrix  of  vectors  centered  around  vector[wn]  */ 
for(j  =  -1  ;  j  <=  1  ;  j++) 

( 

for(k  =  -ix  ;  k  <=  ix  ;  k  +=  ix) 

{ 

iftj  +  k  !=  0) 

( 

iflvector[wn  +  j  +  k][5]  ==  1 ) 


printfT"  %d  %d  %f  %f  %f  %f  \n”,wn,wn  +  j  +  k, 
vector[wn][0],vector[wn  +  j  +  k  ][0],vector[wn][l],vector[wn  +  j  +  k  ][1]); 
V 


) 


sum_mag  +=  vector[wn  +  j  +  k][0]; 
sum_ang  =  vectorfwn  +  j  +  k][  1]; 
sum_vec++; 


) 


) 

iflsum_vec  >=  2) 

( 

vector[wn][0]  =  sum_mag  /  sum_vec; 
vector[wn][l]  =  sum_ang ; 
vector[wn][5]  =  2; 

) 

) 

] 

} 


com  pute_vector() 

{ 
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magnitude  =  vector[wn][0]; 
angle  =  vector[wn][l]; 

xpl  =  vector[wn][2]; 
ypl  =  vector[wn][3]; 

xp2  =  vector[wn][2]  +  magnitude  *  cos(angle); 
yp2  =  vector[wn][3]  +  magnitude  *  sin(angle); 

xp3  =  xpl  +  (magnitude  *  0.75)  *  cos(angle  +  h_angle); 
yp3  =  ypl  +  (magnitude  *  0.75)  *  sin(angle  +  h_angle); 

xp4  =  xpl  +  (magnitude  *  0.75)  *  cos(angle  -  h_angle); 
yp4  =  ypl  +  (magnitude  *  0.75)  *  sin(angle  -  h_angle); 


) 


screen_arrow() 

( 

/* 

i  -  cmap  +  magnitude  *  mratio; 
color(i); 

*/ 

I*  draw  arrow  body  V 

move2i(xpl,ypl); 

draw2i(xp2,yp2); 

/*  draw  arrow  head  */ 
f* 

draw2i(xp3,yp3); 

move2i(xp2,yp2); 

draw2i(xp4,yp4); 

*/ 

) 


paper_arrow() 

( 


/* 

pen_number  =  pratio  *  vector[wn][0]; 
fprintfthp,”SP%d  ;\n”),pen  number); 
*/ 

/*  draw  arrow  */ 

fprintflhp,”  PU%f,%f;  \n”,xpl,ypl); 
fprintRhp,”  PD%f,%f;  \n",xp2,yp2); 

/*  draw  arrow  head  */ 

I* 

fprintflhp,”  PD%f,%f;  \n”,xp3,yp3); 
fprintflhp,”  PU%f,%f;  \n”,xp2,yp2); 
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fprintflhp,”  PD%f,%f;  \n”,xp4,yp4); 
fprintflhp,”  \n”); 


make_cmap() 

( 

/******  make  color  map  ***/ 
k  =  0; 

1  =  256  *  6  /  ncolors; 

/*  black  to  red  V 
red  =  0; 
green  =  0; 
blue  =  0; 

for(  j  =  0;  j  <  255  ;  j  =  j  +  1) 

{ 

red  =  j; 

mapcolor(cmap  +  k, red, green, blue); 
k++; 

) 

/*  red  to  purple  */ 
red  =  255; 
green  =  0; 
blue  =  0; 

for(  j  =  0;j  <  255  ;  j  =  j  +  1) 

{ 

blue  =  j; 

mapcolor(cmap  +  k, red, green, blue); 
k++; 

1 

/*  purple  to  blue  */ 
red  =  255; 
green  =  0; 
blue  =  255; 

for(  j  =  0;  j  <  255  ;j  =j  +  1) 

{ 

red  =  255 -j; 

mapcolor(cmap  +  k, red, green, blue); 
k++; 

) 

/*  blue  to  aqua  */ 
red  =  0; 
green  =  0; 
blue  =  255; 

for(  j  =  0;j  <  255  ;j  =  j  +  1) 

{ 

green  =  j; 

mapcolor(cmap  +  k,red,green,blue); 
k++; 

} 
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/*  aqua  to  green  */ 
red  =  0; 
green  =  255; 
blue  =  255; 

for(  j  =  0;  j  <  255  ;  j  =  j  +  1) 

{ 

blue  =  255  -  j; 

mapcolortcmap  +  k,red,green,blue); 
k++; 

) 

!*  green  to  yellow  */ 
red  =  0; 
green  =  255; 
blue  =  0; 

fort  j  =  0;  j  <  255  ;  j  =  j  +  1) 

( 

red  =  j; 

mapcolortcmap  +  k, red, green, blue); 
k++; 

} 

printfT  Cmap  is  %d  and  cmap  +  k  is  %d  \n”,cmap,cmap  +  k); 
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DISTRIBUTION  LIST 


NUWC-NPT  TM  922051 


Internal 

Codes:  102  (K.  Lima) 

0251 
0261 
0262  (2) 

02244 

22 

38 

3891 

81 

8214  (J.  Castano) 

83 
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